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ABSTRACT 


Adaptive  mesh  refinement  (AMR)  of  hyperbolic  systems  allows  us  to  refine  the  spatial 
grid  of  an  initial  value  problem  (IVP),  in  order  to  obtain  better  accuracy  and  improved 
efficiency  of  the  numerical  method  being  used.  However,  the  solutions  obtained  are  still 
limited  to  the  local  Courant-Friedrichs-Lewy  (CFL)  time-step  restrictions  of  the  smallest 
element  within  the  spatial  domain.  Therefore,  we  look  to  construct  a  multi-rate  time- 
integration  scheme  capable  of  solving  an  IVP  within  each  spatial  sub-domain  that  is 
congruent  with  that  sub-domain’s  respective  time-step  size. 

The  primary  objective  for  this  research  is  to  construct  a  multi-rate  method  for  use 
with  AMR.  In  this  thesis  we  will  focus  on  constructing  a  2nd  order,  multi-rate  partitioned 
Runge-Kutta  (MPRK2)  scheme,  such  that  the  non-unifonn  mesh  is  constructed  with  the 
coarse  and  fine  elements  at  a  two-to-one  ratio.  We  will  use  general  2nd  and  4th  order  finite 
differences  (FD)  methods  for  non-uniform  grids  to  discretize  the  spatial  derivative,  and 
then  use  this  model  to  compare  the  MPRK2  time-integrator  against  three  explicit,  2nd 
order,  single-rate  time-integrators:  Adams-Bashforth  2  (ABo),  Backward  Differentiation 
Formula  2  (BDF2),  and  Runge-Kutta  2  (RK2). 
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I.  INTRODUCTION 


There  is  an  ever  increasing  need  to  develop  numerical  methods  capable  of  solving 
large-scale,  real-world  problems.  In  fact,  many  of  these  problems  arise  in  the  physical, 
biological  and  engineering  sciences,  such  that  we  must  use  numerical  models,  because 
the  problems  cannot  be  solved  analytically.  In  other  words,  we  seek  to  find  numerical 
methods  to  approximate  the  solution  to  a  problem,  which  is  often  modeled  as  a  partial 
differential  equation  (PDE),  or  system  of  PDEs. 

For  example,  the  Department  of  Applied  Mathematics  of  the  Naval  Postgraduate 
School  in  cooperation  with  the  Naval  Research  Laboratory,  both  located  in  Monterey, 
California  are  continually  refining  a  modeling  environment  for  solving  the  three 
dimensional  (3D)  compressible  Euler  and/or  Navier-Stokes  equations;  this  is  a  Non¬ 
hydrostatic  Unified  Model  of  the  Atmosphere,  known  as  NUMA.  The  NUMA  model  is 
constructed  using  an  element-based  Galerkin  framework,  which  allows  either  continuous 
(CG)  or  discontinuous  (DG)  high-order  Galerkin  methods,  and  can  be  used  as  either  a 
mesoscale  or  global  model,  which  exhibits  excellent  scaling  properties.  Additionally,  the 
model  is  extremely  flexible  and  allows  the  user  to  easily  interchange  new  time- 
integrators,  grid  and  data  structures,  or  new  basis  functions,  thereby  ensuring  that  the 
model  is  algorithmically  flexible. 

There  is  a  significant  interest  in  using  the  discontinuous  Galerkin  method  (DG) 
for  solving  fluid  dynamics  problems;  studies  by  [1]  and  [2],  have  shown  that  the  DG 
method  is  a  good  choice  for  the  construction  of  future  non-hydrostatic  numerical  weather 
prediction  models,  as  it  combines  high-order  accuracy  of  the  solution  with  geometric 
flexibility  of  non-confonning  (unstructured)  grids.  In  order  to  increase  the  scale 
resolution  capabilities  of  DG,  as  well  as  to  take  better  advantage  of  available  computing 
power,  the  use  of  adaptive  mesh  refinement  (AMR)  for  a  quadrilateral  non-conforming 
grid,  is  currently  being  investigated. 
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In  Figures  1  and  2,  we  see  a  snapshot  of  a  rising  thermal  bubble  being  solved  on 
both  an  adaptive  and  unifonn  mesh,  using  the  quadrilateral  based  DG  method.  Figure  1 
shows  the  solution  at  time  t  =  0  sec.,  and  Figure  2  is  at  time  t  =  1000  sec. 


adaptive!;!  {{puniform 
mesh  mesh 


Figure  1.  Rising  Thermal  Bubble  (Initial  Time  =  0  sec.) 


Figure  2.  Rising  Thermal  Bubble  (Final  Time  =  1000  sec.) 
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To  motivate  the  usefulness  of  the  multi-rate  approach  let  us  first  describe  the  AMR 
strategy.  In  the  figures  above,  we  note  that  the  meshes  were  refined  every  second 
throughout  a  1000  second  simulation.  For  example,  using  a  time-step  of  At  =  0.005 
would  refine  the  mesh  at  every  200  time-steps,  in  order  to  maintain  the  desired  resolution 
within  the  areas  of  the  grid  where  the  temperature  variation  exists.  The  refined  elements 
were  chosen  at  the  respective  time-steps  by  looking  at  the  value  of  the  temperature  within 
each  element  and  then  determining  whether  that  value  was  below  or  exceeded  a 
predetennined  threshold. 

The  total  runtime  of  the  fully  resolved  solution  (Figure  1,  right  panel)  is 
approximately  287  seconds  compared  to  the  AMR  solution  with  a  runtime  of  88  seconds. 
This  is  a  speed-up  of  nearly  3.26.  Even  though  the  AMR  solutions  require  less  runtime, 
the  overall  numerical  method  is  still  restricted  to  the  smallest  time-step,  which  is 
detennined  by  the  smallest  grid  size  within  the  spatial  domain.  Therefore,  both  the 
uniform  and  adaptive  mesh  solutions  require  the  same  time-step  size.  Although  we  can 
reduce  the  computational  cost  using  AMR,  we  would  like  to  take  advantage  of  the  non- 
conforming  grid  by  simultaneously  using  larger  time-steps  within  the  coarse  regions  of 
the  spatial  domain,  and  smaller  time-steps  within  the  fine  regions.  This  problem  of  using 
a  time-step  size,  commensurate  with  the  element  size,  is  the  motivation  behind  my 
research. 

From  the  literature  we  find  that  the  current  approaches  for  tackling  time- 
integration  include:  fully  explicit,  fully  implicit,  or  a  combination  of  the  two,  known  as 
implicit/explicit  (IMEX)  methods.  Currently,  NUMA  uses  the  IMEX  time-integration 
method,  which  circumvents  the  fast  and  slow  propagating  waves  by  splitting  up  the 
physical  process  of  the  problem  being  solved,  such  that  the  fast  moving  components  are 
handled  implicitly,  while  the  slow  moving  components  are  handled  explicitly.  Although 
the  implicit  solutions  have  no  restriction  on  how  large  of  a  time-step  one  may  take,  the 
overall  numerical  method  is  still  restricted  to  the  time-step  size  detennined  by  the  explicit 
solution’s  spatial  step-size.  Therefore,  we  look  to  develop  a  multi -rate,  time-integration 
method  that  will  allow  multiple  time-steps  to  be  used  simultaneously  within  the  various 
elements  of  the  spatial  domain. 
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A.  OBJECTIVES  OF  RESEARCH 

In  this  thesis,  we  will  focus  on  the  construction  and  development  of  an  explicit, 
second-order,  strong  stability  preserving  (SSP)  Runge-Kutta  (RK),  multi-rate  method, 
which  we  will  use  to  solve  the  first-order,  ID,  advection  equation  with  constant  wave 
speed.  This  multi-rate  method  will  then  be  compared  to  its  equivalent  single-rate  SSP 
RK2  method,  using  both  accuracy  and  efficiency  as  the  two  primary  metrics. 

The  emphasis  of  the  thesis  is  on  the  time-integration  aspects  of  numerical 
modeling,  and  for  this  reason  we  restrict  our  spatial  discretization  to  simple  finite 
difference  (FD)  methods,  albeit  of  high  order,  and  generalized  for  non-uniform  grid 
spacing.  We  will  show  that  the  multi-rate  approach  will  retain  its  formal  order  of 
accuracy,  while  increasing  the  efficiency  of  the  numerical  model.  We  will  show  this 
using  complexity  analysis  (i.e.,  operations  count)  and  measuring  total  wall-clock  time  for 
the  single-rate  versus  multi-rate  approaches. 

This  thesis  is  organized  into  three  main  chapters,  not  including  the  introduction 
and  summary  chapters.  In  Chapter  II  we  will  look  at  how  to  represent  the  spatial 
derivative  of  our  initial  value  problem  (IVP),  by  constructing  a  range  of  finite  difference 
stencils  of  various  orders  of  accuracy.  We  will  then  discuss  how  to  choose  which  stencil 
is  best  applicable  by  looking  at  the  accuracy,  stability  and  convergence  of  each  spatial 
stencil,  in  combination  with  a  forward  Euler  scheme  in  time.  In  Chapter  III,  we  will  shift 
our  focus  to  looking  at  three  different  time-integration  methods  within  the  multi-step  and 
multi-stage  families.  Specifically,  we  will  look  at  the  Adams-Bashforth  (AB),  Backward 
Differentiation  Fonnula  (BDF),  and  Runge-Kutta  (RK)  methods  (each  explicit  in  time) 
for  representing  our  temporal  derivative.  We  will  solve  our  IVP  using  each  time- 
integrator  in  conjunction  with  the  FD  spatial  stencils  developed  in  Chapter  II.  Finally,  in 
Chapter  IV  we  introduce  the  idea  of  a  non-confonning  grid  and  derive  an  explicit, 
second-order,  SSP  RK  multi-rate  method.  Here,  we  use  the  single-rate  results  from 
Chapter  III  in  order  to  compare  against  our  multi-rate  method.  The  thesis  will  be 
concluded  by  a  short  summary  of  findings  and  a  discussion  on  future  research. 


4 


B. 


APPLICATION 


Throughout  this  thesis,  the  governing  equation  we  will  use  as  our  test  case,  is  the 
hyperbolic  PDE 


c\2u  =  0  , 


(1.1) 


which  is  commonly  referred  to  as  the  wave  equation,  with  c  being  the  displacement  from 
rest  (or  wave  speed).  We  chose  to  use  Equation  (1.1)  because  we  can  easily  compute  the 
analytical  solution  using  the  method  of  characteristics,  which  can  only  be  used  for 
hyperbolic  problems  possessing  the  right  number  of  characteristic  families.  Furthennore, 
the  principal  types  of  problems  solved  using  NUMA  are  hyperbolic  PDEs. 

In  order  to  simplify  our  problem,  we  will  assume  that  (1.1)  is  one  dimensional, 
such  that 


2v~7  2 

c  V  u 


d2U  2  S2U 
- T~c  - V 

dt2  dx1 


0  . 


Additionally,  we  can  rewrite  equation  (1.2)  as 


( d  d  \ 

- b  C - 

\dt  dx  J 


'  d  cH 
ydt  dxy 


u  =  0  , 


(1.2) 


(1.3) 


since  the  mixed  derivative  terms,  ±c - ,  cancel  each  other  out.  If  we  allow 

dtdx 


q  = 


r  d  d^ 
Kdt  dx , 


u 


(1.4) 


then  after  substituting  equation  (1.4)  into  equation  (1.3)  we  find  the  first  order,  ID,  wave 
equation,  also  known  as  the  advection  equation  to  be  the  following: 


®L<A=o. 

dt  dx 


(1.5) 


Equation  ( 1 .5)  represents  a  right  moving  wave  for  constant  wave  speed,  such  that  c>  0  . 
Likewise,  we  could  have  also  chosen  to  simplify  equation  (1.3)  as  a  left  moving  wave; 
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however,  we  will  now  consider  Equation  (1.5)  as  the  PDE  we  wish  to  solve  numerically. 
Therefore,  the  IVP  we  wish  to  solve  is  equation  (1.5),  along  with  the  proper  initial  and 
boundary  conditions,  such  that  the  problem  being  solved  is  well-posed. 

For  our  purposes,  and  ease  of  computation,  we  impose  the  following  periodic 
boundary  conditions  for  xe[-l,+l]: 

q(~  1,0  =  ?(U)  (1-6) 

^(-U)  =  ^(U)  ,  (1.7) 

ox  ox 

where  these  boundary  conditions  will  be  helpful  as  we  move  from  our  various  single-rate 
methods  to  the  construction  of  a  multi-rate  partitioned  Runge-Kutta  method,  which  we 
will  see  in  Chapter  IV.  Lastly,  we  know  from  the  literature,  that  the  analytical  solution  of 
our  linear  and  homogeneous  IVP,  Equation  (1.5),  can  be  represented  by  d’  Alembert’s 
solution  where  we  know  that  the  right  moving  wave  can  be  defined  as  q(x,  t)  =  f(x-  ct )  . 
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II.  FINITE  DIFFERENCES 


In  this  chapter,  we  will  look  at  the  method  of  finite  differences  (FD)  using  our 
PDE  from  Chapter  I.  We  will  solve  the  initial-value  problem  (IVP)  using  the  FD 
approach,  where  we  will  change  the  domain  of  the  continuous  problem  into  a  discretized 
fonn.  This  means  that  the  dependent  variable(s)  will  exist  only  at  the  discrete  points  of 
the  grid  generated  by  our  finite  difference  stencil  [3].  Although  we  know  how  to  solve  the 
original  continuous  PDE  analytically,  we  will  show  the  benefit  in  using  finite  differences 
to  estimate  the  solution,  and  how  this  method  can  provide  an  excellent  approximation 
depending  on  the  order  of  accuracy  desired. 

The  reason  for  choosing  a  PDE,  in  which  we  know  the  solution,  is  to  ensure  that 
our  numerical  method  approximates  the  original  PDE  well.  Once  we  know  that  our 
numerical  method  is  consistent  and  accurate,  we  can  then  apply  the  method  to  future 
PDEs  or  systems  of  PDEs  where  we  may  not  necessarily  know  how  to  compute  the 
analytical  solution,  or  finding  the  exact  solution  is  practically  impossible.  In  these  cases, 
we  must  instead  rely  on  numerical  methods  that  utilize  computer  algorithms  to 
approximate  the  solutions  [4]. 

It  is  also  important  to  note  that  in  solving  an  IVP  numerically  with  FD,  we  must 
replace  all  continuous  derivatives  with  finite  difference  approximations;  however,  we  are 
not  required  to  use  the  same  stencil  or  difference  quotient  for  each  derivative  within  the 
problem.  Instead,  we  are  free  to  choose  multiple  finite  difference  schemes  for  each 
derivative  within  the  PDE.  We  will  see  this  later  in  the  chapter;  however,  let’s  first  take  a 
look  at  how  to  build  different  finite  difference  schemes  for  approximating  these 
derivatives. 

A.  SPATIAL  FINITE  DIFFERENCE  APPROXIMATIONS 

Since  the  linear,  first  order,  1-D,  wave  equation  has  only  two  first  order  partial 
derivatives,  one  in  space,  dq/dx,  and  one  in  time,  dq/dt ,  we  only  need  to  choose  two  FD 
stencils.  However,  we  will  explore  other  options  for  time-integration  methods  in  Chapter 
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Ill,  where  we  will  see  that  higher-order  FD  methods  are  not  as  practical  in  time.  For  now, 
we  will  focus  primarily  on  the  spatial  derivative  dq/dx. 

When  approximating  a  derivative  using  FD,  we  must  make  some  basic  choices 
that  involve  how  we  generate  the  grid,  which  will  serve  as  the  discretized  domain  for 
each  solution  in  space  and  time  where  we  want  to  compute  the  solution  of  the  modified 
IVP  [5].  Clearly,  the  finer  the  spatial  grid  spacing,  the  more  solutions  we  are  required  to 
compute  at  each  time  step,  and  the  better  our  approximation  will  be  to  the  original 
continuous  problem;  however,  we  must  consider  the  cost-benefit  analysis  of  higher 
accuracy  versus  computational  cost.  It  is  important  to  note  that  a  FD  stencil  with  higher 
accuracy  is  not  only  a  function  of  the  grid  spacing,  but  also  the  number  of  grid  points 
used  to  approximate  the  derivative. 

Figure  3  shows  a  general  grid  using  unifonn  spacing  in  both  space  and  time.  This 
grid  is  used  to  approximate  the  derivatives  for  the  solution  at  each  grid  point  (xi,tn), 

using  the  grid  point  (xi,tn)  and  its  neighboring  points.  Therefore,  if  we  think  of  the 
solution  to  the  PDE  as  q(x(t),t) ,  then  q"  =  q(xi,tn) . 
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Figure  3.  Unifonn  Grid  in  Space  and  Time 
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Finite  difference  methods  are  simple  to  construct  and  are  beneficial  when  the 
geometry  of  the  IVP  is  regular.  Therefore,  for  our  1-D  wave  equation  we  will  currently 
assume  we  have  a  uniform  grid  with  periodic  boundary  conditions  for  our  analysis. 

B.  FIRST  ORDER  APPROXIMATION 


Using  the  grid  points  in  Figure  3,  we  are  able  to  construct  a  first  order  backward 
FD  approximation  to  the  continuous  spatial  derivative  at  (x.,tn),  using 


q’U  =q(xi-Ax,tn) 

such  that 


dq,  _  g,  -  q ,-_i 
ox  Ax 


(2.1) 


We  begin  with  the  backward-difference  fonnula  instead  of  the  forward-difference 
since  it  is  known  in  advance  that  Equation  (2.1)  is  unconditionally  stable  for  the  IVP 
(stability  analysis  will  be  a  major  topic  throughout  this  thesis  and  will  be  discussed  later 
in  this  chapter,  and  future  chapters,  for  all  numerical  methods  used).  Furthennore,  this  is 
an  obvious  choice  for  the  derivative,  given  that  if  we  take  the  limit  of  the  right  hand  side 
as  Ax  — »  0 ,  we  have 


fix 


li- 1 


Ax— >0 


Ax 


+  0(Ax)  =  lim 

Ax^O 


q(xi>tn)~  q(Xj  -Ax,tn) 
Ax 


+  0{  Ax) 


which  is  the  very  definition  of  the  first  order  derivative  of  q  with  respect  to  x  at  grid 
point (xi,tn) .  Therefore,  we  have  built  a  first  order  approximation  to  the  IVP  such  that  the 

solution  q ”  is  first  order  accurate  (O(Ax))  in  space  for  a  “sufficiently”  small  but  finite 

Ax  [3].  In  fact,  using  Taylor  series  expansion  (TSE)  of  the  backward-difference  formula 
above  will  provide  a  more  rigorous  proof  that  Equation  (2.1)  is  indeed  first  order 
accurate.  This  can  be  seen  in  Appendix  A,  where  we  show  that  the  continuous  derivative 
is  equivalent  to  the  numerical  approximation  plus  the  truncation  error  (T.E.),  such  that  the 
T.E.  =  0( Ax)  . 
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Note  that  the  forward-difference  fonnula 


^  n  n  n 

=  ~9t  +  0(Ax)  (2.2) 

ox  Ax 


is  also  a  first  order  accurate  approximation  to 


The  big  O  notation  used  above  has  a  very  accurate  mathematical  meaning,  such 
that  when  we  replace  T.E.  with  0( Ax) ,  this  represents  the  truncation  error  as  being 

bounded  by  a  positive  real  constant  K,  multiplied  by  the  absolute  value  of  |Ax| ,  where 

At  is  reasonably  small,  such  that  we  can  write  \T.E.\  <  K  |  Ar|  as  At  — >  0 .  In  other  words, 

the  truncation  error  is  on  the  order  of  At  raised  to  the  highest  power  found  within  all  the 
terms  in  the  truncation  error  [3].  We  can  think  of  ( O )  notation  as  the  order  of  accuracy  for 
the  numerical  method,  such  that  if  we  look  at  all  of  the  terms  in  the  TSE,  and  keep  only 
the  tenn  with  the  largest  growth  rate,  then  O(Ax)  means  we  are  keeping  all  terms  of  At  , 

and  throwing  away  all  other  tenns  (At2)  and  higher.  It  is  important  to  note  that  O(Ax) 
does  not  necessarily  inform  us  on  how  large  the  truncation  error  is,  yet  it  does  provide  us 
an  indication  on  how  the  numerical  method  performs  as  At  gets  closer  and  closer  to 
zero.  This  is  accomplished  by  refining  our  grid  spacing. 

For  example,  if  our  numerical  method  is  on  the  order  of  0(Ax2 ) ,  then  as  we 
decrease  At  by  half,  the  estimated  error  of  the  method  will  decrease  by  a  factor  of  four, 
such  that 


O 


f— T 

v  2  J 


=  0 


f 

v 


At2  ^ 


Results  for  this  can  be  found  in  Chapter  III,  where  we  will  look  at  various 
numerical  time-integration  methods  for  solving  our  IVP.  For  more  information  on  ( O ) 
notation,  refer  to  [6]. 
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c. 


HIGHER  ORDER  APPROXIMATIONS 


There  are  many  ways  we  can  develop  FD  approximations  for  dq /dx .  The  nonnal 
convention  for  building  these  representations  begins  by  choosing  the  specific  grid  points 
we  want  to  use  to  approximate  the  derivative,  and  then  expanding  each  point  about  the 
grid  point  (xi,tn) .  The  most  common  first  order  derivative  representations  can  be  found  in 
Table  1,  where  we  only  use  two  to  three  grid  points  to  represent  these  derivatives  [3]. 


Derivative  Finite-Difference  Representation  Equation 

dq'; 

dx 

q‘  q‘-'  +  O(Ax)  (2.3) 

Ax 

qM~qi  +  0(Ax)  (2.4) 

Ax 

+0(Ax2)  (2.5) 

2Ax 

-3q:+4q;+1-q;+  2 

-t-c/^zvx  )  U'Oj 

2Ax 

3<li  ~ 4q,~'  +  q‘-2  +  0(Ax2)  (2.7) 

2Ax 

Table  1.  Difference  Approximations  Using  Two  To  Three  Points 


Equations  (2.3)  and  (2.4)  are  the  backward  and  forward  FD  approximations 
respectively,  whereas  Equation  (2.5)  is  the  centered-difference  approximation  which  can 
be  computed  by  subtracting  the  TSE  of  q(xj-Ax,t")  from  the  TSE  of  q(xt  +  Ax,tn) 

about  the  point  q” .  This  gives  us  a  second  order  accurate  difference  fonnula  for  the  first 

derivative  using  the  points  qnM  and  c/"  , .  Appendix  B  outlines  how  we  can  construct  a  4th 

order,  centered  difference  stencil  for  the  first  derivative,  which  will  be  used  in 
conjunction  with  the  time-integration  methods  developed  in  Chapters  III  and  IV  to 
analyze  single  and  multi-rate  methods  in  both  unifonn  and  non-uniform  grids.  This 
stencil  is  4th  order,  such  that 
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(2.8) 


dgl  q>U-  SC.  +  SC-g” 

fix  12  Ax 


+  0(Ax4)  . 


It  is  easy  to  see  how  we  can  construct  any  FD  stencil  for  first,  second  or  even 
higher  order  derivatives,  as  well  as  mixed  partial  derivatives,  if  desired,  by  using  the 
same  approach  outlined  in  Appendix  B.  However,  the  PDE  we  will  use  only  requires  first 
order  derivatives;  therefore,  we  will  not  need  to  construct  any  more  FD  stencils  than  what 
has  already  been  shown. 


D.  DIFFERENCE  REPRESENTATION  OF  A  PDE 


From  Chapter  I  we  saw  that  the  IVP  in  question  is  the  1-D,  first-order  wave 
equation,  with  constant  wave  speed,  c,  such  that 


dq  _  dq 
dt  dx  ’ 


c  >  0. 


(2.9) 


Using  either  equations  (2.3)  and  (2.4),  or  the  grid  points  found  in  Figure  2,  we  can 
rewrite  our  PDE  in  a  first  order  accurate  discretized  form  for  both  space  and  time.  We 
will  use  forward  Euler  to  represent  the  temporal  derivative,  and  Equation  (2.3),  first  order 
backward-difference  formula,  for  the  spatial  derivative,  such  that 


n  n+ 1  n 

oqj  _  A  -q, 

dt  At 


and 


Ax 


Substituting  these  back  into  Equation  (2.9),  we  then  have  the  following 
approximation  to  our  PDE: 


■q, 


At 


=  -c 


q>  -q,-i 

Ax 


(2.10) 


where  we  can  rewrite  (2.10)  such  that 


q""  =  qni  -v(qni~  qh )  ,  where 


cr  =  c- 


A  t 
Ax 


(2.11) 
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We  will  see  later  on  in  the  chapter,  that  the  value  of  a  (Courant  number)  will 
play  an  important  role  in  determining  where  our  numerical  method  is  stable,  as  it 
measures  how  fast  information  flows  across  the  grid  points. 

Equation  (2. 1 1)  is  the  explicit  forward  Euler  (first-order  upwind)  representation  of 
the  PDE  (2.9)  using  a  first  order  FD  stencil  in  space.  The  method  is  considered  explicit 
since  we  are  solving  the  equation  for  only  one  unknown  value,  qr‘"  ,  which  is  the  solution 

at  the  next  time  step.  As  previously  stated,  the  spatial  and  temporal  FD  schemes  do  not 
have  to  be  the  same;  therefore,  if  we  fix  the  time  stencil  to  be  forward  Euler,  we  have  an 
infinite  number  of  options  for  representing  the  spatial  derivative.  We  could  choose  to  use 
a  higher  order  method  in  space;  however,  the  overall  accuracy  of  the  combined  numerical 
method  will  only  be  as  good  as  the  weakest  link.  This  means  that  if  we  use  a  first  order 
method  in  time  and  a  fourth  order  method  in  space,  the  overall  method  is  still  only  first 
order  accurate.  Below  are  a  few  more  examples  for  the  explicit  forward  Euler  method: 


n+ 1 


<h  -Vi 
At 


=  -c 


g/+i~Wi 

2Ax 


(2.12) 


n+ 1 


<h  ~<li 
At 


=  -c 


3<y;  :  •  q: 

2Ax 


(2.13) 


Both  methods  above  are  second  order  accurate  in  space  such  that  the  T.E.  for  each 
methods  is  0(At,Ax2) .  This  is  important,  given  that  we  can  show  that  the  accuracy  in 
space  and  time  are  independent  of  each  other.  However,  there  are  many  factors  we  must 
consider  when  detennining  which  numerical  approximation  method  to  use  when  solving 
a  particular  problem.  Of  these,  accuracy,  consistency,  stability  and  efficiency  must  be 
analyzed.  Although  we  want  a  numerical  method  that  is  as  accurate  as  possible, 
consistency  and  stability  are  of  much  more  importance,  such  that  if  the  method  is  both 
“consistent  and  stable,  then  it  will  converge  to  the  correct  solution’’  [7].  However,  if  we 
have  two  methods  that  are  both  stable,  consistent,  and  have  the  same  order  of  accuracy, 
then  we  might  want  to  look  at  which  method  is  most  efficient  (i.e.,  lower  computational 
cost).  It  is  important  to  note  that  regardless  of  how  accurate  or  efficient  the  method  is,  if 
we  do  not  have  both  stability  and  consistency,  then  our  numerical  solution  will  be  of  no 
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value  [8].  Therefore,  we  must  first  ensure  the  method  is  stable  and  consistent,  and  then 
we  can  look  for  how  to  increase  the  accuracy  or  efficiency  of  the  method,  depending  on 
what  is  of  most  importance  to  us. 


1.  Consistency 


We  have  looked  at  the  accuracy  of  a  single  FD  stencil;  however  we  have  not 
shown  how  to  determine  the  accuracy  of  the  entire  numerical  method.  Using  Equation 
(2.10),  we  can  show  the  accuracy  of  the  method  by  first  applying  TSE  about  the  point  q " 
and  then  substituting  into  the  numerical  method. 


qf  =q(xi,t  +  At)  =  q"  +  At  +  +  0(At3) 

1  '  '  dt  dt2  2! 


q‘ ,  =  q(x,  -  Ax,  f  )  =  q*  -  Ax +  8  q‘  (Ax) —  6>(Ax3 ) 

M  ‘  '  Sx  Sx2  2! 

Substituting  these  expansions  into  Equation  (2.10)  gives  us  the  following: 


At 


=  -c 


f  n  in  dq"  .  d2q"  (Ax)2  ...  ^ 

q"~  Ax+  -0(Ax  ) 

^  OX  OX  L !  j 

Ax 


where  the  above  equation  reduces  to: 


dt  dt2  2! 


+  0(Ar)  =  -c 


" dq ;  d2q';  Ax 
,  dx  dx2  2! 


+  0(Ax2) 

) 


^  +  c^  +  O(At,Ax)  =  0  .  (2.14) 

dt  dx 

Thus,  we  have  proven  that  method  (2.10)  is  first  order  accurate  in  both  space  and  time. 

Now,  a  numerical  method  is  said  to  be  consistent  if  the  difference  between  its 
approximation  to  the  PDE  and  the  exact  solution,  i.e.,  truncation  error,  goes  to  zero  in  the 
limit,  as  the  grid  is  refined.  Therefore,  if  we  allow  both  the  temporal  grid  spacing 
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parameter,  At ,  and  spatial  grid  spacing  parameter,  Ax ,  in  Equation  (2.14)  to  both  go  to 
zero  in  the  limit,  such  that  (At,  Ax)  — >  0 ,  then  we  can  clearly  see  we  will  recover  the 
original  continuous  PDE,  Equation  (2.9),  thereby  showing  the  numerical  method,  (2.10), 
is  consistent.  This  same  procedure  is  valid  for  any  numerical  method,  where  more 
information  on  consistency  and  stability  for  a  generalized  discretization  can  be  found  in 
most  numerical  textbooks. 

2.  Stability  and  Convergence 

When  deciding  whether  or  not  a  FD  stencil  is  appropriate  to  use  in  solving  an 
IVP,  we  have  looked  at  the  properties  of  accuracy  and  consistency;  however,  one  of  the 
most  critical  properties  to  analyze,  is  the  stability  of  the  method.  In  fact,  not  only 
stability,  but  also  the  rate  at  which  the  numerical  solutions  for  the  method  being  used 
converge  to  the  true  solution  of  the  PDE. 

One  way  to  measure  the  difference  between  the  true  solution,  which  we  will 
define  as  </>” ,  and  our  numerical  solution,  q'‘ ,  is  to  take  the  nonn  of  the  difference 
between  the  two  solutions  at  each  grid  point  [9].  Although  there  are  many  fonns  for 
computing  a  nonn,  we  will  use  the  Euclidean,  or  L2 ,  norm,  which  is  defined  as 

9, 

Therefore,  we  want  to  see  that  in  the  limit,  as  (At,  Ax)  — >  0 ,  the  value  of  the  norm 
between  the  numerical  solution  and  the  true  solution  converges  on  the  order  of 
0(At' ',  Ax')  ,  where  r  and  5  represent  the  order  of  accuracy,  such  that 

\V-q';\2=°(Af,Axs). 

Note  that  this  numerical  method  satisfies  the  Lax  Equivalence  Theorem,  “which 
states  that  if  a  FD  scheme  is  linear,  stable  and  accurate  of  order  (r,s),  then  it  is 
convergent  of  order  ( r,s )”  [9].  In  other  words,  this  theorem  establishes  the  fact  that  if  a 
unique  solution  exists  for  our  IVP,  and  the  solution  depends  continuously  on  both  the 


f  N 
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initial  data  and  boundary  conditions,  then  the  problem  is  considered  to  be  well-posed; 
such  that,  if  we  have  a  consistent  and  stable  method,  then  the  numerical  method  will 
indeed  converge  to  the  true  solution  of  the  original  PDE  [8]. 

The  above  theorem  is  extremely  important  and  forces  us  to  not  simply  develop  a 
FD  scheme  that  is  accurate  or  consistent.  According  to  [9],  consistency  is  not  enough  to 
assure  the  convergence  of  a  numerical  method.  The  method  must  also  be  stable.  There 
are  a  great  number  of  consistent  finite-difference  methods  that  are  utterly  useless  because 
they  are  unstable.  Therefore,  if  we  return  to  the  definition  of  stability,  we  must  ensure 
that  the  numerical  solution,  q " ,  does  not  grow  without  bound,  such  that  for  every  time 
step  nAt ,  less  than  the  final  time  T,  we  can  find  a  constant  CT  that  satisfies 
||<7"||  <  CT  |<7°| ,  for  all  values  of  At  and  Ax  ,  that  are  reasonably  small  [9]. 

3.  Von  Neumann  Stability  Analysis 

We  have  just  seen  the  importance  of  stability;  however,  we  would  now  like  to 
extend  this  idea  to  our  particular  problem  at  hand,  such  that  using  Equation  (2.10)  or 
(2.11),  we  can  find  for  what  values  of  a  our  method  will  be  stable.  Any  value  that 
satisfies  the  stability  conditions  will  then  fall  into  what  we  call  the  stability  region. 

The  most  widely  known  and  used  procedure  for  analyzing  stability  is  the  Von 
Neumann  Method,  such  that  Von  Neumann’s  stability  analysis  looks  at  the  discretized 
solution  for  a  given  time  step,  and  represents  this  solution  as  a  finite  Fourier  series,  where 
the  overall  stability  of  the  numerical  method  can  be  detennined  by  analyzing  each 
element  or  Fourier  mode  of  the  series  [9].  In  this  manner,  we  can  ensure  that  the  overall 
stability  of  our  method  will  indeed  be  stable  if  every  Fourier  mode  is  stable. 

Below  is  the  Fourier  series  representation  for  the  numerical  solution  q'j , 

«;  =  £ 

1--N 
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where  q\“]  represents  the  amplitudes  of  the  wave,  i  =  V-l ,  k ,  are  the  wave  numbers  and 
k,  Ax  represents  the  phase  angle,  sometimes  denoted  (ft, ,  which  detennines  the  low  from 
high  frequency  waves,  such  that  a  phase  angle  of  zero  indicates  a  large  amplitude  wave 
and  a  phase  angle  of  n  indicates  a  small  amplitude  wave  [8]. 

It  would  be  extremely  tedious  to  evaluate  all  Fourier  modes;  therefore,  if  we 
instead  look  at  just  one,  where  we  ensure  that  this  mode  has  an  amplitude  q"\  that  does 
not  grow  for  all  time,  and  is  bounded  such  that  |g(n)|  <  1,  then  we  can  determine  where 

our  region  of  stability  lies.  For  more  detailed  infonnation  on  Von  Neumann’s  Method, 
look  to  [3],  [8],  and  [10], 

If  we  now  concentrate  on  the  single  Fourier  mode 

n  ~{n)  i]<t> 

q,  =  q  eJ  , 

then  we  can  substitute  this  mode  into  Equation  (2.1 1),  with  the  following  statements: 

q"  =  Gneij* ,  q"+l  =  Gn+1eu * ,  q)_x  =  G  V0_1) * , 

such  that  we  have 

G"+lej*  =  Gneyt  -a (Gne‘Ml  - Gne‘u~m) 

where  Gn  replaces  q  "  ]  as  the  amplitudes.  If  we  now  simplify  the  above  expression  by 
dividing  through  by  Gn  e1'1’ ,  we  find  that 

G  =  l-o-(l-<f"). 

Using  Euler’s  fonnula,  e"1’  =  cos (f> ± ism (f) ,  offers 

G  =  1  -  cr(l  -  cos  ^)-/cr  sin  (2.15) 

which  we  know  is  the  graph  of  a  circle  centered  at  (1  -  cr,  0)  with  radius  cr .  Note  that  we 
want  to  solve  Equation  (2. 15)  for  cr  values  that  satisfy  |G|  <  1  where 
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(2.16) 


Ax 

From  Equation  (2.15)  we  find  that  the  solution  is  a  <  1  and  is  referred  to  as  the 
Courants-Friedrichs-Lewy  (CFL)  condition,  which  ensures  that  our  solution  for  this 
specific  numerical  method  will  be  stable  within  the  stability  region  defined  by  cr  <  1  [10]. 
This  can  be  seen  in  Figure  4.  This  figure  depicts  the  regions  of  stability  for  Equation 
(2.11),  such  that  stability  is  maintained  for  values  inside  of  the  unit  circle  (solid  blue 
line),  given  a  particular  Courant  number.  Several  values  for  cr  have  been  provided, 
which  demonstrates  that  for  0  <  cr  <  1 ,  the  method  is  stable. 

Another  method  for  analyzing  the  stability  region  is  to  look  at  the  magnitude  of 
the  amplification  factor  |G|,  such  that  for  G  <  1  the  method  is  stable.  Figure  5  shows  a 

contour  plot,  where  each  contour  line  represents  a  value  of  |G|  for  0.5  <cr<  1.15  and 
0  <  (f>  <  In .  Here  we  chose  a  few  values  of  cr  >  1 ,  in  order  to  verify  that  they  indeed 
correspond  to  values  of  |G|  >  1 ,  therefore  proving  that  Equation  (2.11)  is  unstable  for 
Courant  values  greater  than  one. 


Stability  Region  for  Eq  (2.11) 


Re(G) 


Figure  4.  Stability  Region  for  Forward  Euler  w/  1st  Order  Upwind 
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Stability  Region  for  Eq  (2. 11) 


Figure  5.  Contour  Plot  for 


Euler  w/  1st  Order  Upwind 
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III.  TIME-INTEGRATION  METHODS 


In  Chapter  II  we  introduced  finite  difference  (FD)  methods  and  demonstrated  how 
we  can  use  these  methods  to  approximate  spatial  derivatives.  FD  are  based  on  Taylor 
series  expansion  (TSE),  and  are  also  used  to  derive  time-integration  stencils,  which  we 
will  use  to  solve  our  initial  value  problem  (IVP).  In  this  chapter,  we  will  turn  our  focus  to 
two  branches  of  time-integration  methods  known  as  multi-step  and  multi-stage.  We  will 
not  go  into  great  detail  for  either  of  these  methods,  as  they  can  be  found  in  most 
numerical  analysis  texts.  However,  we  will  briefly  introduce  these  methods  for  the 
purpose  of  comparing  each  method  with  single-rate  results  using  our  particular  IVP. 

A.  MULTI-STEP  METHODS 

Multi-step  methods  are  all  methods  that  require  starting  values  from  several 
previous  time  steps.  These  methods  can  be  very  useful;  however,  they  require  prior 
infonnation  from  the  system  being  solved.  Therefore,  they  must  first  be  started  by  using 
an  initial  value,  then  implementing  a  one-step  method  to  receive  two  values,  then  a  two- 
step  method  to  receive  three  values,  and  so  on,  depending  on  the  order  of  the  method.  In 
fact,  Forward  Euler  (explicit  method)  and  Backward  Euler  (implicit  method)  are  first- 
order  multi-step  methods,  although  they  both  contain  only  one-step. 

Other  multi-step  methods  include  the  Trapezoidal  and  Leapfrog  schemes,  which 
are  both  implicit,  second-order  accurate  methods.  In  general,  there  are  two  particular 
families  of  multi-step  methods:  the  Adams  methods  and  the  Backward  Differentiation 
Formulas.  We  will  start  by  looking  at  the  Adams  methods. 

1.  Adams  Methods 

Within  the  Adams  family  of  multi-step  methods,  the  two  most  commonly  used  are 
Adams-Bashforth  of  order  p,  (AB/<),  and  Adams-Moulton  of  order  p,  ( AM/<),  For  our 
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purposes,  we  want  to  only  focus  on  explicit  time  integration  methods;  therefore,  we  will 
only  look  at  A  B />,  as  these  methods  are  explicit  in  time,  whereas  the  AMp  are  all  implicit 
in  time. 

The  general  fonnula  for  the  Adams-Bashforth  method  for  solving  the  IVP, 
q'  =  F(t,q),  q(t0)  =  q0,  t>t0, 

is 

m 

q”'  ,  (3.1) 

j= 1 

such  that  any  AB/>  method  can  be  derived  by  Equation  (3.1),  which  can  be  constructed 
using  either  Taylor  series  expansion  or  Lagrange  interpolating  polynomials.  If  we  turn 

3  1 

our  focus  to  AB2,  where  m  =  2  ,  b,  =  — ,  and  b0  =  —  using  Equation  (3.1),  we  can  see 

2  2 

that  we  achieve  the  following  formula, 

qn+1=qn  +y(3  F{tn,qn)-F(t*-\qn-%  (3.2) 

which  is  a  two-step  method,  and  requires  the  solutions  at  two  previous  times,  q"  and 
q"  1 ,  in  order  to  compute  the  next  solution,  qn"  ,  which  is  at  time  tn+1 .  Therefore,  we  can 
use  any  single-step  method  we  like  to  jumpstart  AB2;  however,  in  order  to  maintain  2nd 
order  accuracy,  the  single-step  method  must  also  be  2nd  order  accurate.  For  a  more 
detailed  description  of  how  Equation  (3.2)  is  constructed  using  TSE,  look  to  Appendix  D. 

Before  using  AB2  with  our  IVP,  it  is  necessary  that  we  know  the  stability  region 
of  this  method;  therefore,  if  we  first  assume  that  there  is  no  error  within  our  spatial 
discretization  method,  then  if  we  let  F(tn  ,qn)  =  Xq" ,  and  F(t"  \qn  1 )  =  Aq"  1 ,  then  when 
we  perfonn  the  Von  Neumann  stability  analysis  for  AB2,  we  have  the  following: 
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where  z  =  A  At 


n+1 

q 


i(n+l)6  _  gine 


|(  3eM -eKn~l)e' 


where  we  let  q"  =  e"'° ,  because  we  seek  to  find  the  curve  for  which  \q\  =  1 ,  and  defines 
the  boundary  of  stability.  Then,  solving  for  z,  we  find: 

2(cos  6  +  i  sin  6)  -  2 

Z  —  - . 

3  -  cos  6  +  i  sin  6 

If  we  now  separate  the  real  from  the  imaginary,  we  find  that: 


Re(z)  = 


-2  cos2  6  +  8  cos  9  +  2sin 0  -  6 
(3  -cos6)2  +  sin2  0 


T  ,  .  -4  sin  9  cos  6  +  8  sin  6 

Im(z)  = - - - - — 

(3  -cos  9)  +sin26 

such  that,  plotting  the  Re(z)  versus  the  Im(z)  provides  the  region  of  stability  for  AB2, 
which  is  shown  in  Figure  6.  From  the  figure,  we  see  that  for  any  value  of  z  =  A  At  chosen 
within  this  region,  the  solutions  to  the  IVP  will  remain  stable.  However,  as  stated 
previously,  this  is  assuming  we  have  zero  error  within  our  spatial  discretization  method. 
This  is  important  as  it  allows  us  to  see  where  each  time-integration  method  is  stable, 
independently  of  the  spatial  scheme. 


23 


AB2  Stability  Region  Assuming  No  Spatial  Error 


Figure  6.  AEF  Stability  Region 


We  have  now  seen  the  general  fonnulas  for  Forward  and  Backward  Euler,  and 
AB2.  In  fact,  using  the  following  multi-step  formula, 

m— 1  m— 1 

<?"*'  =I>/T'  +A  tY.bftt"1  ,r‘)  ,  (3.3) 

7=0  7=-l 

if  we  let  m  =  1 ,  then  we  have  a  one-step  method,  such  that  Equation  (3.3)  becomes, 
qn+l  =  a0qn  +  At  (b_lF(tn+\q',+1)  + b0F(tn  ,qn)) . 

Now,  if  we  let  a0  =  b0  =  1 ,  and  b_x  =  0  ,  we  find  that  the  above  equation  becomes 

qn+l  =qn  +AtF(tn,qn), 

which  is  the  exact  fonnula  for  the  explicit  Euler  Method.  Therefore,  we  see  that  Euler  is 
indeed  a  multi-step  method  for  /u-step  equal  to  one,  and  is  equivalently  an  Adams- 
Bashforth  method  of  order  one  (ABi). 

Now,  using  Equation  (3.3)  and  the  coefficients  from  Table  2,  we  can  formulate  a 
few  of  the  most  commonly  used  multi-step,  time-integration  methods  of  order  p  =  1,2,3 . 
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Coefficients 


Method  (P) 

Euler  (1) 

sis 

Backward  Euler  (1) 

Trapezoidal  * (2) 

Leapfrog  (2) 

BDF2  *  (2) 

AB2  (2) 

AM2  *  (3) 

Table  2.  Multi-step  Methods  of  Order  p  =  1,2,3  *Implicit  Method 

From  Table  2  we  also  notice  that  depending  on  the  value  of  the  b;j  coefficient,  for  j  =  1 , 

then  Equation  (3.3)  will  generate  either  explicit  or  implicit  time-integration  methods.  In 
other  words,  for  all  values  of  b_x  ^  0  ,  then  the  multi-step  method  will  be  implicit  in  time. 

2.  Backward  Differentiation  Formulas 

Another  family  of  multi-step  methods  is  known  as  the  Backward  Differentiation 
Formulas  (BDF)  of  order  p,  such  that  BDF2  represents  a  two-step  method,  which  is 
second  order  accurate  in  time.  It  is  important  to  note  that  the  number  of  steps  in  the 
method  does  not  necessarily  determine  the  order  of  accuracy  of  the  method.  In  general,  it 
is  possible  to  prove  that  the  maximal  order  of  a  convergent  implicit,  in- step  method,  using 
Equation  (3.3),  is  at  most  in  +  2  if  m  is  even,  and  in  +  1  if  m  is  odd;  however,  for  explicit 
schemes,  an  /u-step  method  cannot  attain  an  order  greater  than  in  [11],  This  is  known  as 
the  first  Dahlquist  barrier.  For  example,  Table  2  lists  AM2  as  a  two-step  method; 
however,  this  is  an  implicit  time-integration  method  with  3ld  order  accuracy  in  time. 


a0=b0=  1,  b  ,  =  0 


a0  =  1,  b_x  =1,  b0=  0 


ao  =  h  b  i  =b0=~,  bt=  0 
2 


a0  =  0,  a,  =  1,  b_ j  =  hj  =  0,  b0  =  2 


4  l  2 

ao  ~  —  ’  a\  —  -  ~  ’  b_  j  —  — ,  b0  —  bx  —  0 


3  l 

a0  —  a,  —  b_ .  =  0,  b0  =  — ,  b .  =  — 
2  2 


a0  =  h  ,  ax  =  0,  b_x  =  — ,  b0  =  — ,  bx  =■ 
12  12 


1 

12 
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Although  Equation  (3.3)  is  the  general  equation  for  all  multi-step  methods,  we 
saw  that  there  was  a  more  general  expression  for  the  family  of  AB/>  methods  using 
Equation  (3.1).  This  is  also  true  of  BDF/>,  such  that  we  can  express  every  BDF/>  method 
using  the  following  equation: 

m 

'^Jam_Jqn+1~j  =  AtF(tn+1,qn+1)  ,  (3.4) 

j=o 

Previously,  we  showed  that  ABi  was  equivalent  to  Forward  Euler  (explicit 
method).  In  a  similar  fashion,  we  can  show  that  Backward  Euler  (implicit  method)  is 
equivalent  to  BDFi,  for  m  =  1,  a0  =  ax  =  1 .  Yet,  of  more  interest  to  us,  is  BDF2,  as  this 

method  will  provide  us  with  an  additional  2nd  order  time-integration  method  to  use  for 
our  IVP.  However,  given  that  we  want  to  only  focus  on  explicit  time-integrators,  we  must 
look  at  how  we  can  rewrite  BDF2  so  that  it  is  instead  explicit  in  time. 


First,  let  us  consider  BDF2: 

h+1  d  n  1  n_\  2 At  ,  .n+l  m+1\ 

q  =  - q  --q  +— ^  )  (3-5) 

We  can  see  that  the  implicit  form  of  BDF2  requires  infonnation  from  two  previous  time 
steps,  in  addition  to  the  solution  at  time  tn+l .  In  order  to  solve  for  qn+l ,  we  would  need  to 
solve  a  system  of  equations;  therefore,  in  order  to  remove  the  implicitness  of  the 
equation,  we  can  instead  use  an  approximation  for  F(q"" )  by  the  method  of 
extrapolation.  Since  we  can  solve  the  IVP  for  the  solutions  at  F(q” )  and  F(q“  1 ) ,  then  if 
we  say  that 

A Fn  =F(qn)-F(q"-1)  and  A Fn+l  =  F(qn+1)- F(qn)  , 
and  then  approximate  AF”+1  «  A F" ,  we  will  find  that 

F(qn+1)~ F(q")  «  F(q")~ F(q"~l) 

F(qn+1)*2F(qn)-F(q n  l)  . 
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However,  this  approximation  uses  the  assumption  that  we  can  use  a  linear 
extrapolation.  Therefore,  if  we  approximate  F(qn+l)  by  F(q" )  and  F(q"  1 )  using  TSE 
about  F(qn+1 )  ,  we  find: 

(a)  F(q")  =  F(q"+1)-dF(f+1)  At  +  0(At2) 

ot 

(b)  F(qn~l)  =  F(q"+1)~  2  8F )  At  +  0(At2)  . 

dt 

Now,  if  we  multiply  (a)  by  2  and  subtract  (b),  then  we  get  an  expression  for  F(q"+1) , 
such  that, 

F(qn+l)  =  2F(q")-F(q'1-1)  . 

Although  both  methods  provided  the  same  expression  for  F(q"+1) ,  we  have  found  that 
the  solution  using  the  TSE  approach  is  exact,  whereas,  we  cannot  always  guarantee  the 
linear  extrapolation  method  will  work.  If  we  now  substitute  this  approximation  back  into 
Equation  (3.5),  we  have  an  explicit  formula  for  BDF2,  such  that 

q  =-q  ~-q  +—{2F(t  ,q  )-F(t  ,q  ))  .  (3.6) 

Refer  to  Appendix  D  for  a  more  thorough  derivation  of  Equation  (3.6). 

Let  us  now  look  at  the  Von  Neumann  stability  analysis  for  Equation  (3.6),  where 
we  again  assume  there  is  no  spatial  error,  such  that, 

?“I=V"_V"I+T 

where  we  let  q "  =  e",(> ,  q"  =e',n  l)/y,  and  qn+l  =e,(n+1>0,  such  that  factoring  out  the  term 
e‘nB  provides  the  following  expression: 

ew=---e-w  +—(2-e-ie). 

3  3  3  v  ’ 

Now,  solving  for  z,  we  find  that 
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z  = 


— cos#  — 
v3  3  j 


4^\  7 

+  z  —  sin  9 
3 


4  2  , 

- cost/ 

3  3 


A 


+  i  — sin£? 
3 


Therefore,  separating  the  real  from  the  imaginary,  we  finally  have  the  following: 


Re(z)  = 


Im(z)  = 


-2cos  d  +  6cos9  +  sin  9  -  4 
cos2  6  -4  cos  9  +  4 

-3  sin  9  cos  t?  +  4  sin  <9 
cos2  9  -4cos#  +  4 


such  that,  a  plot  of  the  stability  region  for  BDFt  is  shown  in  Figure  7,  where  we  see  that 
for  any  value  of  z  =  XAt  chosen  within  this  region,  the  solutions  to  the  IVP  will  also 
remain  stable  assuming  no  spatial  error. 


BDF2  Stability  Region  Assuming  No  Spatial  Error 


Figure  7.  BDFt  Stability  Region 
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B.  MULTI-STAGE  METHODS 

Like  multi-step  methods,  multi-stage  methods  have  the  desirable  property  that 
they  can  achieve  very  high-order  accuracy,  while  simultaneously  reducing  the  amount  of 
derivatives  that  need  to  be  computed  at  each  grid  point.  In  order  to  design  a  multi-stage 
method,  let  us  first  consider  the  ordinary  differential  equation  (ODE),  y'  =  F(t,y) ,  such 

that  if  we  integrate  from  time  tn  to  t"+l ,  then  we  can  show  that, 

tn*  i  t "+1 

y(t"+l)-y(tn)=  J  y\t)dt  =  \  F(t,y{t))dt . 

tn  tn 

If  we  now  apply  the  trapezoidal  rule,  we  get 

y(tn+1)-y(tn)  =  ^(F(t\y(0)  +  F(tn+\y(tn+l)))  + 
which  is  the  Trapezoidal  method  in  Table  2,  such  that, 

n+ 1  n  .  /  77»/vA2  n  \  .  rp/^+l  n+l\\ 

y  =y  +—{F{t  ,y  )  +  F{t  ,y  )), 

where  the  global  error  is  0(At2) ;  therefore,  this  method  is  a  second  order,  implicit  multi- 
step  method.  It  is  important  to  note,  that  in  general,  where  we  have  written  At ,  we 
typically  write  h,  such  that  h  is  any  step  size,  and  is  not  restricted  to  only  time.  However, 
for  our  purpose  in  solving  the  IVP,  we  will  use  these  methods  as  time-integrators  where 
we  let  h  =  At . 

Although  we  have  an  expression  for  an  implicit,  multi-step  method,  we  instead 
want  an  explicit,  multi-stage  method.  To  accomplish  this,  we  will  now  use  Forward  Euler 
to  solve  the  ODE  for  v'i+1 ,  such  that  Euler’s  method  will  be  called  a  “predictor,”  and 
labeled  as  v"+1 .  We  then  substitute  yn+l  into  the  Trapezoidal  method,  known  as  the 
“corrector,”  which  gives  us  the  following: 

Euler:  yn+1  =  y"  +  AtF(t"  ,yn)  (predictor) 

Trapezoidal:  v”+1  =y"  +^-^F(t"  ,yn)  + F(tn+l ,  y'!+1))  (corrector) 
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The  above  combination  of  Euler  and  Trapezoidal  methods  is  referred  to  as  Heun’s 
method,  and  is  now  an  explicit,  second  order,  multi-stage  method.  We  may  also  choose  to 
look  at  the  above  method  by  writing  it  the  following  way: 


K  =  A  tF(t",y”) 
k2=A tF(tn+1,y"  +kl)  , 

yn+1  =yn+^{kl+k2) 


(3.7) 


where  the  number  of  stages  is  represented  by  k,  ,  (i  =  2 .  Here,  Heun’s  method, 
Equation  (3.7),  is  known  as  a  two-stage  method,  for  p  =  2  . 

1.  Explicit  Runge-Kutta  Methods 

We  have  just  shown  how  to  construct  a  very  simple,  explicit,  multi-stage  method, 
where  Heun’s  method  is  classified  within  a  much  larger  family  of  multi-stage  methods 
known  as  Runge-Kutta  of  order  M  (RKm).  Here,  Heun’s  method  is  most  commonly 
referred  to  as  RK2. 

In  general,  any  explicit  p-stage  RKm  method  can  be  constructed  in  the  following 

way: 


kx=AtF{t\yn) 


p 


where  any  RKm  method  is  determined  by  the  coefficients  p,  c,  d,  and  b ,  which  are 
typically  displayed  in  a  table  referred  to  as  a  Butcher  tableau,  after  J.C.  Butcher  [10]. 


dT  c 


b 
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p  -  Number  of  stages.  This  is  not  to  be  confused  with  the  order  of  the  method; 
hence,  change  of  subscript  to  M  in  RKM.  For  example,  there  are  fifth  order  RK 
methods  with  six  stages. 

c-  Ap  xp  coefficient  matrix. 

d  -  Row  vector  of  size  p. 

b  -  Row  vector  of  size  p. 

Although  we  will  not  focus  on  the  derivation  of  how  these  parameters  are  chosen  for 
specific  RKm  methods,  the  reader  can  find  more  information  in  “Numerical  Methods  for 
Ordinary  Differential  Equations,”  by  Butcher  [10]  or  “Numerical  Analysis,”  by  Burden 
and  Faires  [12]. 

Below  are  examples  for  RK2,  RK3  and  RK4,  with  their  associated  Butcher 
Tableau: 


RK2 


kl=AtF(t\y”) 

<  k2  =AtF(tn+\yn  +kx) 

y'"=y"+\(K+k2) 


kx  =  AtF(t"  ,yn) 

k2=AtF(tn+^-,y"+^ 

k3=AtF(t"+At,yn+kl  +2k2) 

yn+'=yn+-k  +  1-k2+-L 
l  o  3  o 


0 

0 

0 

1 

1 

0 

1 

“T 

2 

2 

0 

1 

1 

2 

2 

1 

1 

2 

1 

2  1 

6 
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kx  =  AtF(tn,yn) 

(  A  t  k  \ 

k2=AtF  f+  —  ,y"+^ 

{  2  2) 

k,  =  Atc[V +y,y"  +  yj 

k4  =  AtF(tn  +  At,  y"  +  k3 ) 

y"+i  =y"  +-(kx+2k2  +2k3+k4) 

6 

Recall  from  Chapter  II,  if  our  numerical  method  is  of  order  0(hp),  then  as  we  reduce  the 
step  size  by  a  factor  of  two,  the  estimated  error  of  the  method  is  reduced  by  a  factor  of 
2P  .  Therefore,  for  any  /;-th  order  method,  if  we  let  p  =  M  and  h  =  At ,  then  for  RK4, 
which  has  global  truncation  order  of  0(At4) , 


In  general,  both  multi-step  and  multi-stage  methods  have  their  advantages  and 
disadvantages.  Higher  order  multi-step  methods  require  more  use  of  past  values,  whereas 
multi-stage  methods  require  more  calculations  per  step.  For  RK  methods,  the  main 
computational  effort  is  in  the  evaluation  of  the  right  hand  side  (RHS)  function  itself  [10]. 
For  example,  in  the  second-order  RK  methods,  the  local  truncation  error  is  0(Ar), 
where  the  cost  is  two  functional  evaluations  per  step,  given  there  are  two  stages  within 
the  step  [10],  Likewise,  the  fourth-order  RK  methods  require  four  functional  evaluations 
per  step,  as  there  are  four  stages,  and  the  local  truncation  error  is  0(At4).  However,  this 
pattern  does  not  hold  for  all  RK  methods.  In  general,  Table  3  shows  “the  relationship 
between  the  number  of  evaluations  per  step  and  the  order  of  the  local  truncation 
error”  [10]. 


0 

1 

1 

2 

2 

1 
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1 
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Number  of  Evaluations  /  Step  and  Order  of  the  Local  Truncation  Error 


Evaluations  per  step  2 

3  4 

5  <  ft  <  7 

8<  «  <9 

10  <  ft 

Best  possible  local 

truncation  error  0(At2) 

0(At3)  0(At4) 

0{Atn~x) 

0(At"~2) 

0(At"3) 

Table  3. 

Evaluations  per  Step  and  Truncation  Error 

Given  the  result  from  Equation  (3.8),  it  is 

easy  to  see 

why  a  higher  order  RK 

method  might  be  preferred  to  a  lower-order  RK,  if  we  are  concerned  with  accuracy; 
however,  Table  3  shows  that  in  comparing  higher-order  RK  methods  ( n  >5)  with  lower- 
order  RK  methods  ( n  <  4  ),  that  the  lower-order  methods  may  instead  be  preferable  to 
higher-order.  In  other  words,  although  the  higher  order  methods  allow  for  a  larger  time 
step  and  improved  local  truncation  error,  the  overall  computational  cost  increases 
significantly.  Therefore,  depending  on  the  problem  being  solved,  we  may  be  satisfied 
with  using  a  3ld  or  4th  order  RK  method,  which  requires  a  slightly  smaller  time  step, 
compared  to  say,  an  8th  order  method;  however,  the  overall  number  of  functional 
evaluations  per  step  is  reduced  drastically.  There  is  always  a  tradeoff  between  choosing  a 
method  that  has  the  best  possible  local  truncation  error,  yet  minimizes  cost. 

For  our  analysis,  we  will  focus  only  on  explicit,  second-order  methods;  therefore, 
let  us  now  look  at  the  stability  of  RKo,  which  will  be  our  primary  time-integration  method 
for  evaluating  the  IVP  in  Chapter  I.  If  we  rewrite  Equation  (3.7)  in  the  following  way: 

qn+l=qn+^{F{qn)  +  F{qn+AtF{qn)j)  ,  (3.9) 

then,  after  applying  Von  Neumann  stability  analysis  to  Equation  (3.9),  we  find  the 
subsequent  representation  for  z,  such  that, 

z2  +  2z  -  2e'6  +2  =  0.  (3.10) 

Here,  Equation  (3.10)  is  quadratic;  therefore,  we  must  now  solve  for  the  roots  of  the 
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function,  and  then  plot  the  real  versus  imaginary  values  of  z,  for  0  <  6  <  2n  .  Figure  8 
displays  these  results,  where  we  now  see  the  stability  region  for  Equation  (3.9)  assuming 
no  spatial  error. 


RK2  Stability  Region  Assuming  No  Spatial  Error 


Figure  8.  RK2  Stability  Region 

Although  these  plots  are  useful  on  their  own,  it  is  more  beneficial  if  we  overlay 
each  method’s  stability  region,  in  order  that  we  may  compare  AB2,  BDF2,  and  Rfo 
against  each  other.  Figure  9  shows  the  stability  plot  for  each  of  these  methods,  in  addition 
to  first-order  Euler,  where  we  can  easily  see  that  RK2  has  a  significantly  larger  region  of 
stability  in  comparison  to  the  other  three  explicit  methods.  Also,  this  plot  visually 
supports  the  fact  that  RK2  is  stable  for  much  larger  time-step  values  than  the  other  three 
methods. 
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Stability  Regions  Assuming  No  Spatial  Error 


Figure  9.  Stability  Regions  for  AB2,  BDF2,  and  RK2 


C.  SINGLE-RATE  RESULTS  ON  UNIFORM  GRID 

Let  us  now  look  at  the  numerical  results  for  solving  our  IVP  on  a  uniform  grid, 
using  the  multi-step  and  multi-stage  time-integrators  mentioned  in  sections  A  and  B. 
First,  if  we  choose  our  initial  condition  to  be  a  Gaussian  function  of  the  form: 

-( x-bf 

q{x)-ae  2cl  ,  (3.11) 

where  a,  b,  and  c  are  real  positive  constants,  such  that, 
a:  height  of  the  curve’s  peak, 
b:  center  position  of  curve’s  peak, 
c:  width  of  bell  curve, 
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then  we  know  the  general  graph  of  this  function  is  a  symmetric  bell  curve  that  has  tail 
ends,  which  fall  off  to  positive  and  negative  infinity.  Figure  10  shows  three  Gaussian 
functions  with  varying  parameters  to  demonstrate  how  each  parameter  affects  the  shape 
of  the  curve. 


Gaussian  Functions  for  Various  Parameters 


x 


Figure  10.  Gaussian  Plot  for  Various  Parameters 

For  our  analysis  of  the  IVP,  we  will  set  the  domain  to  be  {x  e  M  |  -1  <  x  <  1} . 
Additionally,  we  will  impose  periodic  boundary  conditions,  such  that  the  solutions  at  the 
right  and  left  boundaries  are  equal  to  each  other  as  the  1-D  wave  equation  propagates  to 

the  right  in  time,  for  0  <  t  <  1 ,  and  A t  =  a  ^ ,  where  the  time-step  is  a  function 
detennined  by  the  Courant  number,  wave  speed  and  grid  spacing.  For  the  following 
results,  we  will  set  the  wave  speed  to  be  constant  at  c  =  2 ,  a  will  vary  depending  on  the 
stability  of  the  numerical  method  and  the  grid  spacing,  A* ,  will  be  detennined  by  the 
number  of  grid  points  used  to  evaluate  the  solution  for  a  particular  time-step. 
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1.  Explicit  RK2  Method 


Let  us  first  rewrite  our  IVP,  Equation  (2.9), 


dq  _  dq 
dt  dx 


c  >  0, 


using  Equations  (2.5)  and  (3.7),  where  (2.5)  is  a  second-order  centered  difference 
approximation  for  the  continuous  spatial  derivative,  and  (3.7)  is  the  RK2  time-integration 
method.  Here,  we  know  that  if 


Cs  n  n  n 

F(q" )  =  -c-^-  =  -c  q‘  '  +  0( Ax- ),  for  c  =  2, 


dx  2Ax 

then,  substituting  into  (3.7)  we  find: 


Ax  v  7 


k2  =  A  tF{q.  +kl)  =  A  tF 

At 

Ax 


n  In  n  \  | 


n  A  t  ( 


A  t,  n 


Therefore,  we  now  have  the  following  expression  for  RK2  using  a  second-order  centered 
difference  stencil  in  space: 


+  - 


At 

2Ax' 


-2q, ”+?■,). 


(3.11) 


If  we  fix  our  initial  parameters  to  be  cr  =  0.5  and  Ax  =  0.02 ,  then  when  we 
approximate  Equation  (2.9)  using  Equation  (3.11),  Figure  11  shows  a  graphical 
comparison  between  the  exact  solution  and  the  numerical  FD  solution  at  times  t  =  0.025, 
0.50,  1.0,  and  4.0. 
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Solution  using  RK2  at  time  =  0.25  Solution  using  RK2  at  time  =  0.5 


Figure  11.  Numerical  vs.  Exact  Solution  for  RKo  using  2nd  Order  CFD  with  Ar  =  0.02 


At  time  t  =  0,  the  numerical  solution  (red  curve)  is  set  equal  to  the  initial 
condition  (blue  curve),  and  as  time  increases,  the  initial  condition  (Gaussian  function) 
begins  to  propagate  to  the  right.  From  Figure  1 1  we  can  clearly  see  that  after  one 
revolution,  where  t  =  1,  the  difference  between  the  numerical  solution  and  initial 
conidtion  is  approximately  0.0249,  and  at  t  =  4,  the  difference  is  0.0986;  where  these 
nonn  values  are  computed  using  the  L2  norm,  which  was  defined  in  Chapter  II  as 


f  N  \/2 

Zkf  . 

V  !=1  J 


such  that  the  difference  between  the  numerical  solution,  q" ,  and  exact  solution,  <f>" ,  is 
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{  N  \X 

’  for  n  =  i’2’3’- 

V  i=l  J 


It  is  easy  to  see  from  Figure  11,  that  the  numerical  solution,  defined  by  the 
parameters  above,  quickly  loses  accuracy  as  time  increases,  such  that  both  the  dispersion 
and  dissipation  errors  are  evident.  The  dispersion  error  is  a  result  of  the  numerical 
solution  evolving  slower  in  time  compared  to  the  real  solution,  while  the  dissipation  error 
is  the  difference  in  the  height  of  the  numerical  solution  to  that  of  the  real  solution. 
Therefore,  since  At  is  a  function  of  the  Courant  number,  wave  speed  and  grid  spacing, 
we  have  a  few  options  to  take  in  order  to  determine  if  we  can  achieve  better  results  for 
solving  (2.9)  using  (3.11):  we  can  increase  or  decrease  the  number  of  grid  points  used  to 
evaluate  the  solution  (i.e.,  make  Ax  smaller  /  larger),  modify  the  wave  speed,  choose  a 
different  Courant  number,  or  a  combination  of  the  three,  so  long  as  we  ensure  that  the 
stability  of  (3.11)  is  maintained. 

Let  us  look  at  the  results  for  only  changing  At  .  If  we  change  the  number  of  grid 
points  from  100  to  400,  then  At  =  0.005,  and  our  time-step  now  becomes,  At  =  0.00125  , 
compared  to  the  previous  time-step  value  of  At  =  0.005 ,  for  a  fixed  Courant  number  of 
cr  =  0.5  .  Running  the  same  simulation  again  with  RfC  for  0  <  t  <  4  ,  we  notice  from 
Figure  12  that  the  numerical  solutions  are  more  accurate  when  using  At  =  0.005,  as 
opposed  to  At  =  0.02  for  t  <  2 ;  however,  for  t  >  2 ,  the  solution  not  only  loses  accuracy, 
but  is  unstable  for  the  given  parameters.  Therefore,  we  see  that  although  the  initial 
computations  using  a  finer  grid  provide  better  results  for  t<  2 ,  the  end  result  is  worse. 

In  fact,  what  we  find  is  that  the  numerical  solution  is  unstable  for  RK2  using  a 
fixed  Courant  value  of  cr  =  0.5  and  a  At  =0.005.  Additionally,  the  computational  cost 
also  increased  as  we  computed  values  at  400  points  versus  only  100  for  each  time-step. 
Initially,  this  did  not  make  sense  as  I  assumed  using  a  finer  grid  spacing  and  smaller  time 
step  would  produce  better  results.  However,  what  we  find  is  that  as  we  refine  the  grid 
spacing  smaller  and  smaller,  the  time  step  value  needed  to  ensure  stability  must  also  get 
smaller. 
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RK2  results  for  Ax  =  0.005  .  At  =  0.00125  .  L2  Norm  =  0  0015498,  T  =  1  RK2  results  for  Ax  =  0  005  ,  At  =  0  00125  ,  L2  Norm  =  0  0031038.  T  =  2 


RK2  results  for  Ax  =  0.005  ,  At  =  0.00125  ,  L 2  Norm  =  0.070785,  T  =  3  RK2  results  for  Ax  =  0.005  ,  At  =  0.00125  ,  L2  Norm  =  32.3054,  T  =  4 


Figure  12.  Numerical  vs.  Exact  Solution  for  RK2  using  2nd  Order  CFD  with  Ax  =  0.005 

Therefore,  we  need  to  look  at  the  other  parameters,  namely  a ,  and  find  a  solution 
that  achieves  the  accuracy  we  desire,  while  simultaneously  reducing  the  computational 
cost  and  ensuring  stability.  Furthermore,  the  results  above  show  the  importance  of 
carrying  out  the  simulation  for  increasing  time.  If  we  had  stopped  the  simulation  after 
only  one  or  two  revolutions,  we  would  have  assumed  the  solution  for  RK2  using  the  finer 
grid  for  the  parameters  chosen  was  better,  when  in  fact  we  have  shown  the  opposite. 
Regardless,  neither  solution  is  desirable. 

Instead,  let  us  look  at  the  solution  to  Equation  (3.11)  by  fixing  the  grid  spacing 
and  only  varying  the  Courant  value.  If  we  now  plot  the  numerical  solution  using  the 
Courant  value  and  estimated  error,  we  find  that  as  the  Courant  value  decreases,  the 
estimated  error  of  the  solution  decreases.  These  results  are  shown  in  Figure  13,  and  they 
confirm  Equation  (3.8),  which  shows  that  as  the  time-step  is  reduced  by  a  factor  of  two, 
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the  estimated  error  of  the  method  is  reduced  by  a  factor  of  2P .  Here,  RFC  is  a  second 
order  method;  therefore,  as  At  is  reduced  by  half,  the  estimated  error  should  decrease  by 
a  factor  of  four.  Although  we  are  plotting  against  the  Courant  value,  remember  that 

At 

<j  =  c  — , 

Ax 

therefore,  there  is  a  direct  correlation  between  a  and  At ,  such  that  as  a  gets  smaller, 
At  must  also  get  smaller,  given  a  fixed  wave  speed  and  grid  spacing. 


Figure  13.  Estimated  Error  vs.  Courant  for  RK2  using  2nd  Order  CFD  with  Ax  =  0.02 


From  Figure  13,  we  notice  that  the  estimated  error  begins  to  increase  for  Courant 
values  of  <7  <  1 .  What  we  find,  is  that  this  is  not  due  to  the  time-integrator,  but  rather  the 
spatial  error  introduced  by  our  approximation  of  the  spatial  derivative.  Here,  we  used 
Equation  (2.5),  which  is  the  second-order,  centered  difference  stencil.  For  our  purposes, 
we  are  not  as  concerned  with  the  spatial  error,  and  only  want  to  focus  on  the  errors 
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produced  from  the  various  time-integration  methods.  Therefore,  if  we  use  a  spatial 
discretization  stencil  with  higher-order  than  the  order  of  our  time-integrators,  we  will 
hopefully  be  able  to  better  analyze  the  solutions  for  decreasing  time-steps  or  smaller 
Courant  numbers.  For  that  reason,  let  us  now  modify  the  numerical  method  so  that  we 
instead  compute  the  spatial  derivative  using  Equation  (2.8), 


fix  12  Ax 


+  0(Ax4) . 


which  is  the  4th  order,  centered  difference  stencil  constructed  in  Appendix  B. 

If  we  again  look  at  the  estimated  error  versus  Courant  value  plot  using  RK2  in 
time  and  Equation  (2.8)  in  space,  then  we  find  the  following  results,  which  are  shown  in 
Figure  14. 


Figure  14.  Estimated  Error  vs.  Courant  for  RK2  using  4th  Order  CFD  with  Ax  =  0.02 
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From  Figure  14,  we  now  notice  that  our  estimated  error  has  improved  drastically  for 
decreasing  values  of  cr ;  however,  our  solution  becomes  increasingly  unstable  for  values 
of  cr  approximately  larger  than  0.75.  Figure  14  clearly  demonstrates  the  usefulness  of 
higher-order  spatial  discretization  methods  when  analyzing  lower-order  time  integrators, 
such  that  we  are  better  able  to  see  the  convergence  rates  of  the  time-integrator,  as  well  as 
the  estimated  error.  For  example,  both  Figures  13  and  14  show  the  estimated  error  for  the 
entire  numerical  method.  Therefore,  we  notice  that  the  spatial  error  dominates  in  Figure 
13,  so  that  we  are  unable  to  see  the  temporal  errors;  whereas  in  Figure  14,  the  spatial 
error  is  not  readily  seen  until  we  reach  a  Courant  value  of  approximately  0.1,  thereby 
allowing  us  to  more  accurately  analyze  the  temporal  error  for  decreasing  cr . 

If  we  now  vary  both  the  grid  spacing  and  Courant  values,  the  results  are  shown  in 
Figure  15.  Here  we  find,  that  as  the  grid  size  decreases,  the  range  of  cr  values  must  also 
decrease  in  order  to  maintain  stability.  This  plot  also  demonstrates  the  trade-off  between 
accuracy  and  efficiency,  such  that  as  Ax  decreases,  the  more  accurate  the  RK2  solutions 
are  to  the  exact  solution;  however,  the  number  of  computations  per  step  must  increase. 


Figure  15.  Error  vs.  Courant  for  Various  Ax 
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Figure  15  is  also  useful  in  helping  to  better  explain  the  previous  RK2  results,  which  were 
shown  in  Figures  1 1  and  12.  Since  we  were  using  a  fixed  Courant  value  of  0.5  and  only 
varying  Ax  ,  we  notice  that  the  blue  curve  supports  the  results  shown  in  Figure  1 1,  as  we 
found  that  although  the  method  gradually  incurred  both  dissipation  and  dispersion  errors, 
it  was  still  stable.  However,  for  the  same  Courant  value  of  0.5,  we  notice  that  the  green 
curve  represents  the  results  shown  in  Figure  12,  such  that  Figure  15  supports  the  fact  that 
for  Ax  =  0.005,  the  numerical  method  is  unstable. 

Similar  to  Figure  15,  it  is  also  useful  to  plot  the  time-step,  At,  versus  the 
estimated  error  of  the  numerical  method.  In  Figure  16,  we  notice  the  areas  along  each 
curve,  where  the  slopes  are  approximately  equal  to  2.  It  is  within  these  regions  that  for  a 
particular  Courant  value  and  spatial  grid  size,  that  the  overall  numerical  method  RK2  is 
stable  and  2nd  order  accurate.  The  analytical  proof  that  RK2  is  second-order  accurate, 
regardless  of  the  spatial  discretization  method,  can  be  found  in  Appendix  C. 


Figure  16.  Error  vs.  Time-step  for  Various  Ax 
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Figure  16  is  also  helpful,  as  it  clearly  shows  the  relationship  between  the  time- 
step  size  and  estimated  error,  as  opposed  to  the  Courant  value.  It  is  easily  seen  that  for  the 
given  spatial  grid  sizes,  the  estimated  error  is  reduced  as  Ax:  decreases;  however,  it  is  at 
the  cost  of  a  much  smaller  time-step,  which  requires  more  computational  time. 

In  this  chapter,  we  presented  three  different  time-integration  schemes:  two  multi- 
step  (ABt  and  BDF2),  and  one  multi-stage  (RK2).  Therefore,  let  us  also  present  the  results 
solving  our  IVP  for  all  three  methods. 


Figure  17.  Error  vs.  Time-step  for  RK2,  AB2,  and  BDF2 


From  Figure  17,  we  notice  that  each  of  the  three  time-integrators  converges  to  the  actual 
solution  at  the  expected  rate  of  approximately  two  as  each  method  is  0(At2 ,Ax4)  .  We 
also  notice  that  RK2  is  stable  for  larger  time-step  sizes,  which  supports  the  stability 
analysis  plots  shown  in  Figure  9.  Therefore,  given  the  above  results,  we  will  use  RK2  as 
the  time-integration  scheme  of  choice  for  the  remainder  of  this  thesis  and  in  developing 
our  multi-rate  method  in  Chapter  IV. 
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IV.  MULTI-RATE  METHODS 


In  this  chapter,  we  focus  on  the  development  of  a  second-order  multi-rate 
partitioned  RfG  method  (MPRK2),  which  uses  a  series  of  convex  combinations  of  Euler 
steps  [13].  To  construct  this  method,  we  begin  by  introducing  a  non-uniform  grid,  where 
we  must  generalize  the  spatial  derivative  stencils  developed  in  the  previous  chapters. 

A.  NON-UNIFORM  GRIDS 

In  Chapters  II  and  III,  we  solved  our  IVP  using  second  and  fourth-order  accurate 
FD  spatial  discretization  methods  that  were  analyzed  on  a  uniform  grid.  However,  as  we 
begin  to  construct  our  multi-rate  time-integration  method,  we  must  now  introduce  a  non- 
uniform  mesh  in  both  space  and  time.  For  now,  we  will  simplify  the  problem  to  only  non¬ 
uniformity  in  space. 


t 


Ac  Ac  Ax  Ac  Ac 

2  2  2 

Figure  18.  Non-unifonn  Grid  in  Spatial  Domain 
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Using  Figure  18,  we  notice  that  this  grid  is  uniform  in  time  for  constant  At ,  yet 
non-uniform  in  space.  Here,  we  let  every  grid  point  to  the  left  of  xt ,  for  all  time,  be  a 
distance  of  Av  apart,  while  all  points  to  the  right  of  xi  are  a  distance  of  Ax/ 2  apart.  It  is 

important  to  note,  that  in  general,  we  can  arbitrarily  vary  the  distance  between  any  two 
points  within  the  non-uniform  grid,  and  that  Figure  18  is  only  one  graphical  example  of 
how  we  may  choose  to  define  our  grid  space.  Furthennore,  this  figure  demonstrates  that 
if  we  now  want  to  build  a  discrete  approximation  to  the  spatial  derivative  of  our  IVP, 
then  we  will  not  be  able  to  use  the  spatial  stencils  developed  in  Chapter  II  and  Appendix 
B,  as  they  all  rely  on  a  single  Ar  . 

For  example,  Equations  (2.5)  and  (2.8), 


'An  n 

_qM-q, 


i- 1 


dx 


2Ax 


+  0(  Ax2) 


dq"  ..  W 2-8gM+8g,”+i-g; 

dx  12Av 


+2 


+  0(  Ax4). 


are  the  second  and  fourth-order  centered  difference  stencils,  respectively,  which  we  used 
to  approximate  the  spatial  derivative  of  our  IVP  in  Chapter  III.  Let  us  now  look  at  how  to 
construct  a  general  expression  for  these  stencils  using  the  grid  points  at  time,  t" , 
referenced  in  Figure  18.  We  will  begin  with  the  second-order  centered  difference  stencil, 
such  that  we  want  to  use  the  solutions  at  q/+v  q/  and  q/ , ,  to  approximate  the  first-order 

derivative  of  q‘‘  with  respect  to  x.  Therefore,  we  now  write  the  spatial  discretization  as  a 
weighted  sum  of  these  solutions,  where 


d£ 

dx 


=  a.q'L  +  P,ql + aw 


(4.1) 


We  will  begin  with  using  TSE  about  q/ ,  such  that, 
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n  n  ( 

qM=q<  +-rM*,+i 
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0[(x,.-x;_i)3_  . 


If  we  now  substitute  these  expansions  into  Equation  (4.1),  then  we  achieve  the  following 
expression: 
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dx 
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7 

Since  we  want  an  expression  for  the  second-order,  centered  difference  stencil,  we 
require  the  following  equalities: 
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=  0 
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which  can  also  be  written  as  the  following  matrix  problem, 
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where  we  let  Ap  =  (xM  -  xi )  and  Am  =  (xt  -xM)  .  After  solving  for  the  coefficients  using 
Cramer’s  rule,  we  find  the  following  expressions: 


a;  = 


A„ 


A  A  2  +  A  A  2  ’ 

pm  m  p 


Pt  = 


-(a  2-a  2) 

V  m _ p  ) 

A  A  2  +  A  A  2  ’ 

p  m  m  p 


-A 


r,  = 


A  A  2  +  A  A  2 

pm  m  p 


Now,  if  we  substitute  the  values  for  the  coefficients  an  (in  and  ;/  into  Equation  (4.1),  we 
find  a  general  expression  for  the  spatial  derivative,  such  that 


M. 

dx 


A  A 

p  n 


+  A  A 


-  + 


o(a,a.)  . 


(4.2) 


which  depends  on  the  solution  at  three  points  with  two  distinct  A  parameters. 

In  order  to  verify  that  Equation  (4.2)  is  the  correct  general  fonnula  for  a  second- 
order,  centered  difference  stencil,  then  if  we  assume  that  both  A  parameters  are  equal 
(i.e.,  unifonn  grid),  we  can  easily  show  that  when  Ap  =Ax  =  Am,  then  Equation  (4.2) 
simplifies  to  Equation  (2.5), 


dq” 


dx 


n  n 


g<+i-g/-i 

2Ax 


+  0(Ax2) , 


thereby  proving  that  Equation  (4.2)  is  a  second-order  accurate  stencil. 

Likewise,  if  we  want  to  construct  a  general,  fourth-order  centered  difference 
stencil,  then  following  the  same  methodology  as  shown  for  developing  Equation  (4.2), 
we  now  find  that  the  spatial  derivative  becomes 


M 

dx 


=  a,q',\-2  +  P,q'L  +  r,q + pqU  +  n.q, 


-2  9 


(4.3) 


where  after  using  TSE  about  q" ,  we  achieve  the  following  matrix  problem: 
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such  that, 

A„2=(*/+2-*/)>  AP1=(XM-X,) 

Am2=(Xi-Xi-2)’  A,„l=(Xi-Xi-l)  ■ 

Note  that  the  general  expression  for  the  matrix  above  is 
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where  k  is  the  order  of  the  method,  such  that  k  is  even. 

After  solving  the  matrix  problem,  we  find  the  following  expressions  for  the 
coefficients: 
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(4.3.1) 
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(4.3.3) 

(4.3.4) 
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Therefore,  if  we  substitute  (4.3.1)  -  (4.3.5)  into  (4.3),  we  achieve  the  fourth-order, 
centered  difference  stencil  we  desired,  such  that 


dqi 

dx 


-  ai<li+2  +  PAm  +  Yfl!  +  %"-l  +  1li(li-2+0(AmAm2A pA pi)  ' 


(4.4) 


Notice,  if  Apl  =  Aml,  Ap2  =  Am2,  and  2Apl  =  A  2,  then  we  can  easily  show  that  Equation 
(4.4)  is  equivalent  to  Equation  (2.8),  such  that 


n  n  o  n  ,  o  n  n 

dq,  _q<-2-H-i+H+x-qi+2 


dx 


12  Ax 


+  0(Ax4). 


Let  us  return  to  our  IVP  ,  where  we  now  solve  for  solutions  using  the  non-uniform 
grid  in  Figure  18.  If  we  use  RKt  as  our  time-integrator  and  Equation  (4.2)  as  the  spatial 
derivative,  then  we  achieve  the  following  results,  which  are  shown  in  Figure  19: 
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Figure  19.  Non-uniform  RKo  Results:  L2  Error  Norms  vs.  Courant  Value  using  2nd  Order 

Centered  Finite  Difference  Method 

From  Figure  18,  we  notice  that  our  domain  is  partitioned  into  two  sub-domains,  such  that 
the  left  sub-domain  represents  the  course  grid  and  the  right  sub-domain  corresponds  to 
the  fine  grid.  Within  Figure  19,  we  demonstrate  that  if  we  fix  the  fine  grid,  Axfme  =  0.02 , 

and  vary  Axcoarse  from  0.02  to  0.08 ,  then  as  the  width  of  Axcoaraf,  increases,  the  overall 
estimated  error  of  the  numerical  method  increases.  We  also  notice  that  when  we  set 
Ax coarSe  =  fine  >  we  return  to  a  uniform  grid,  such  that  the  errors  are  consistent  with  the 

results  shown  in  Chapter  III,  specifically,  Figure  13.  This  makes  sense,  as  we  should  not 
expect  to  achieve  better  results  using  the  non-uniform  grid,  compared  to  the  uniform  grid, 
provided  the  highest  resolution,  Ax/hte ,  within  both  grids  is  equivalent. 

It  is  also  important  to  analyze  the  L2  error  norms  of  both  the  uniform  and  non- 

uniform  results,  such  that  we  should  expect  the  same  convergence  rates.  For  example,  we 
have  shown  for  the  uniform  grid,  that  if  our  time-integration  method  is  on  the  order  of 
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0(At2),  then  as  we  reduce  the  step  size  by  a  factor  of  two,  the  estimated  error  of  the 
method  is  reduced  by  a  factor  of  2P  ,  such  that  for  RK2,  where  p  =  2  ,  we  know 


Va^ 

2 

=  0 

fA'2l 

[Uj 

l  4  J 

However,  when  solving  on  the  non-uniform  grid,  the  time-step  is  now  detennined  by  the 
highest  spatial  resolution,  in  order  to  ensure  stability  of  the  overall  method;  therefore,  At 
is  a  function  of  Axfme .  Figure  20  displays  the  L2  error  norms  for  our  uniform  and  non- 

uniform  grid  results  using  RK2  and  Equation  (4.2).  From  this  figure,  it  is  easily  seen  that 
for  the  unifonn  grids,  as  Ax  decreases  by  half,  the  estimated  error  of  the  method 
decreases  by  a  factor  of  four;  likewise,  for  the  non-uniform  grids,  as  Axfme  decreases  by 

half,  where  Axcoarse  =  2Axfme ,  the  overall  method  also  decreases  by  a  factor  of  four.  Figure 

20  also  verifies  that  the  non-uniform  grid  results  achieve  larger  errors  than  the  unifonn 
grid,  when  both  methods  use  the  same  high  resolution.  For  example,  the  non-uniform  grid 
using  Axfme  =  0.02  and  Axcoarse  =  0.04 ,  achieves  better  results  than  the  uniform  grid  with 

Ax  =  0.04 ;  however,  is  worse  than  the  unifonn  grid  with  Ax  =  0.02 . 


Figure  20.  Unifonn  vs.  Non-uniform  RK2  Results 
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It  is  also  helpful  to  analyze  the  estimated  errors  of  the  unifonn  vs.  non-uniform 
RK2  results,  using  the  degrees  of  freedom  of  the  spatial  mesh.  Therefore,  Figure  21 
demonstrates  that  as  the  degrees  of  freedom  of  the  numerical  method  increase,  the 
estimated  error  decreases,  such  that  the  uniform  grid  will  achieve  better  results  than  the 
non-uniform  grid,  given  an  equivalent  number  of  points.  This  is  the  price  one  must  pay 
for  using  the  geometric  flexibility  of  non-uniform  grids,  which  must  be  used  to  yield 
more  efficient  solutions. 


Figure  2 1 .  Estimated  Error  vs.  Number  of  Points  for  Unifonn  and  Non-unifonn  RK2 

Results 


Although  Figures  20  and  21  are  helpful  in  comparing  the  uniform  and  non- 
uniform  meshes,  we  notice  that  Figure  20  primarily  shows  the  spatial  error  for  the  overall 
numerical  method,  which  is  on  the  order  of  0(At2  ,Ax2) .  Therefore,  in  order  to  better 
analyze  the  temporal  errors,  we  will  now  use  RK2  in  combination  with  Equation  (4.4),  the 
general  fourth-order  centered  difference  stencil,  to  solve  our  IVP,  such  that  the  overall 
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numerical  method  is  on  the  order  of  0(At2 ,  Ax4) ,  when  we  set  Axcoarse  =  Axfme .  Figure  22 

displays  these  results,  where  we  are  now  able  to  more  accurately  analyze  the  temporal 
errors  for  RKo. 


RK2  Results  Comparing  Uniform  and  Non-Uniform  Grids  for  4th  Order  Centered  Difference  Stencil 


Figure  22.  Unifonn  vs.  Non-Uniform  RK2  Results  using  a  4th  Order  CFD  Stencil 


From  Figure  22,  we  clearly  see  that  the  unifonn  results  are  equivalent  to  those  found  in 
Figure  16,  and  as  we  expected,  introducing  the  fourth  order  centered  difference  stencil  in 
space,  has  allowed  us  to  view  the  temporal  errors  more  effectively.  We  also  notice  from 
Figure  22,  that  our  single-rate  RK2  method  performs  as  we  expected.  For  example,  just  as 
we  found  in  Figure  20,  the  results  above  show  that  for  a  given  non-unifonn  grid,  the 
relative  errors  are  greater  than  the  uniform  grid  errors,  when  both  methods  use  the  same 
high  resolution.  Since  we  have  validated  that  our  general  second  and  fourth-order  spatial 
schemes,  in  conjunction  with  our  single-rate  RK2  method  are  both  consistent  and 
accurate,  let  us  now  look  at  how  to  develop  our  multi-rate,  partitioned,  RK2  time- 
integrator. 
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B. 


MULTI-RATE  GRID  DEVELOPMENT 


In  the  previous  section  we  looked  at  the  single-rate  RK2  results  on  a  non-uniform 
grid,  such  that  the  non-uniformity  only  existed  within  the  spatial  domain,  as  was  shown 
in  Figure  18.  Let  us  now  extend  this  non-confonning  grid  to  both  space  and  time,  where 
Figure  23  displays  a  graphical  representation  for  this  non-uniformity. 


/ 


i-2  /-I  i  z+1  z+2  /+3 


Ax  Ax  Ax  Ax  Ax 

2  2  2 

Figure  23.  Non-confonning  Grid  in  both  Space  and  Time 


From  Figure  23,  we  now  notice  that  the  non-unifonnity  within  the  temporal  domain  will 

present  a  problem  in  constructing  our  finite  difference  representation  of  the  continuous 

spatial  derivative  for  specific  points.  Therefore,  if  we  imagine  that  every  grid  point  to  the 

left  of  point  “i”  in  Figure  23  is  of  dimension  At  ,  and  that  every  grid  point  to  the  right  of 

point  “i”  has  dimension  Ar/2 ,  then  it  can  be  easily  seen,  that  if  we  want  to  represent  the 

spatial  derivative  of  our  IVP  using  the  general  fourth-order  centered  difference  stencil 

developed  in  Chapter  IV,  that  we  will  have  no  issues  in  constructing  this  stencil  for  grid 

points  to  the  left  of  “i,  ”  or  right  of  /  +  1 ,  assuming  that  the  grid  extends  infinitely  in  both 
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directions,  thereby  excluding  any  boundary  conditions.  However,  we  notice  in  Figure  24 
that  as  time  increases,  we  lack  knowledge  of  information  at  every  half  time-step,  that  is 
necessary  for  building  the  fourth-order  centered  FD  stencil  for  grid  points  “i”  and  i  + 1 . 
These  locations  are  indicated  by  red  markers  in  Figure  24. 
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Ax  Ax  Ax  Ax  Ax 
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Figure  24.  Lack  of  Information  Needed  to  Build  a  4th  Order  CFD  Stencil 


Let  us  now  introduce  some  tenninology,  such  that  we  will  commonly  refer  to  the 
coarse  region  within  our  spatial  domain  as  the  “slow  region,”  and  the  fine  region  within 
the  same  domain  as  the  “fast  region,”  where  the  grid  point  located  on  the  boundary 
between  these  two  regions  is  an  “interface  point.”  Let  us  also  assume  from  this  point 
forward  that  when  we  refer  to  a  FD  stencil,  that  the  reference  stencil  is  the  fourth-order 
centered  FD  stencil  found  in  Equation  (4.4)  with  weighted  coefficients  defined  by 
Equations  (4.3.1)  -  (4.3.5).  From  Figure  24,  we  notice  that  in  order  to  represent  the 
spatial  derivative  at  grid  points  “i”  and  i  + 1 ,  that  we  must  make  an  approximation  for  the 
values  at  each  half  time-step  where  we  do  not  have  any  information. 
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Using  the  information  depicted  in  Figure  24,  if  we  choose  to  let  the  values  at  grid 
points  x i  t,tn+ 2  and  x i2,t  "+2  be  equal  to  xt_vt”  and  x i  2,t"  respectively,  then  it  can  be 
easily  shown  using  TSE,  that  this  approach  will  reduce  the  order  of  accuracy  of  the 
overall  method.  Likewise,  if  we  choose  to  let  these  same  points  be  equal  to  x  it,tn+]  and 

x  i-2 ’  t  ”+' ,  then  this  naive  approximation  will  also  reduce  the  overall  order  of  accuracy. 
However,  if  we  know  how  to  compute  the  information  for  grid  points  i  —  1  and  i  —  2  at 
time  levels  t "  and  t ,!+I ,  then  a  better  approach  would  be  to  take  an  average  of  the  two 
values  and  use  this  information  as  placeholders  for  approximating  our  derivative,  such 
that  we  can  easily  show  using  TSE,  that  this  averaging  approach  will  maintain  the  desired 
level  of  accuracy,  such  that  our  FD  scheme  will  remain  0( Ax4)  . 

For  our  particular  IVP,  we  are  concerned  with  the  boundary  conditions,  given  that 
we  enforce  periodicity  at  these  boundaries.  Therefore,  in  order  to  simplify  the  problems 
associated  with  a  non-confonning  grid  in  both  space  and  time  at  these  boundaries,  let  us 
define  a  new  grid  that  consists  of  three  sub-domains  (coarse  /  fine  /  coarse).  Figure  25 
shows  a  graphical  representation  of  this  grid,  where  the  first  and  third  sub-domains  (slow 
regions)  consist  of  same-sized  elements,  and  the  second  sub-domain  (fast  region)  consists 
of  elements  half  the  size  of  the  slow  regions. 


1st  Domain  2nd  Domain  3rd  Domain 

(slow  region  1)  (fast  region)  (slow  region  2) 
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Figure  25.  Non-conforming  Grid  with  Three  Sub-domains 


This  grid  is  different  from  the  two  sub-domain  grid  in  Figure  24,  in  that  we  now  have  two 
interfaces  which  we  will  define  as  a  slow/fast  interface  and  a  fast/slow  interface,  and  can 
be  seen  in  Figure  26.  In  addition,  we  also  have  what  we  will  define  as  “buffer  regions.” 
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These  buffer  regions  exist  at  each  interface,  such  that  the  size  of  a  particular  buffer  region 
is  dependent  upon  the  grid  points  required  to  construct  the  FD  stencil.  First,  we  must 
construct  slow  and  fast  grids,  which  will  be  used  to  compute  the  slow  and  fast  numerical 
solutions  at  a  given  time  level.  From  Figure  26  we  notice  that  the  buffer  regions  are 
defined  as  the  intersections  of  the  slow  and  fast  grids.  Note,  it  is  important  to  distinguish 
between  the  overall  coordinate  grid  and  the  three  sub-grids,  as  both  slow  grids  in  domains 
1  and  2,  each  share  grid  points  with  the  fast  grid.  Notice  also  that  the  fast  and  slow  grids 
identified  in  Figure  26  are  not  the  same  as  the  fast  and  slow  sub-domains  shown  in 
Figure  25. 

Slow  /Fast  Fast  /Slow 

Interface  Interface 

\  / 

i  i  i — — — i  i  nun  •••  min  i  i — — — i  i  i  - 

1 - r - - r - A - r - 1 

Ax  Ax/2  Ax- 

Interfaces 


Slow  /  Fast  Fast  /  Slow 

Interface  Interface 


Buffer  Regions 


Figure  26.  Interfaces  and  Buffer  Regions  for  Non-Conforming  Grid 

Now  that  we  have  identified  the  difference  between  the  slow  and  fast  sub- 
domains,  as  well  as  the  slow  and  fast  grids,  we  must  also  define  where  our  actual 
numerical  solution  will  exist,  such  that  in  order  to  compute  the  entire  solution  at  a  given 
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time  level,  we  require  both  the  fast  and  slow  solutions,  where  the  individual  solution  at 
any  given  grid  point  “i,  ”  is  located  in  one  of  three  solution  regions:  purely  fast,  purely 
slow,  or  at  an  interface.  However,  the  fast  solution  will  take  care  of  the  interface  solution. 
Let  us  refer  to  Figure  27  to  see  what  we  mean.  From  this  figure,  we  notice  that  the  fast 
solution  is  computed  on  the  grid  points  within  the  fast  sub-domain  (including  interface 
points)  defined  in  Figures  25  and  26,  whereas  the  slow  solution  is  computed  using  the 
grid  points  within  the  two  slow  sub-domains  minus  the  interface  points. 
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Figure  27.  Computing  a  Fast  or  Slow  Solution  in  the  Buffer  Region 


For  example,  if  we  want  to  compute  the  FD  fast  solution  at  the  first  interface  point 
located  within  the  slow/fast  buffer  region,  then  we  require  infonnation  from  the 
neighboring  grid  points,  such  that  if  the  interface  point  is  grid  point  “i,  ”  then  we  also 
need  the  values  at  grid  points  i-2,  i-1,  z+1,  and  i  +  2,  such  that  the  values  at 
i  -  2,  and  /  - 1  will  come  from  the  slow  solution  and  the  values  at  i  + 1,  and  i  +  2  will 
come  from  the  fast  solution.  However,  as  we  saw  in  Figure  24,  the  slow  solutions  are 
computed  on  the  slow  domain,  such  that  no  information  is  stored  at  every  half  time-step. 
Therefore,  let  us  now  look  at  how  to  construct  the  multi-rate  partitioned  RK2  (MPRKo) 
time-integrator  that  will  help  us  overcome  this  issue. 
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c. 


MULTI-RATE  PARTITIONED  RK2  METHOD  (MPRK2) 


As  we  begin  to  construct  our  MPRK2  scheme,  we  want  to  ensure  that  all  of  the 
previous  properties,  such  as  consistency,  stability  and  accuracy,  as  discussed  in  Chapter 
III  also  hold  for  this  time-integration  technique.  According  to  [13],  if  we  consider  the 
second  order,  Runge-Kutta  method,  RK2,  defined  by  the  Butcher  tableau  as 


kl=AtF(tn,qn ) 

0 

0 

0 

RK2:< 

k2  =  A  tF  [tn  +  At,  q"  +k^ 

1 

1 

0 

q"+1  =q"  +^(kl+k2) 

1 

1 

2 

2 

then  using  the  following  compact  notation  for  explicit  Euler  steps, 

s( At, q)  =  q(t )  +  A tF(t, q)  ,  (4.6) 

as  defined  in  [13],  we  can  rewrite  the  above  RK2  method  as  a  linear  combination  of  Euler 
steps,  such  that  we  can  guarantee  that  the  above  method  will  strongly  preserve  stability. 
This  is  known  as  a  strong  stability  preserving  (SSP)  method,  such  that  SSP  time- 
integrators  are  methods  that  ensure  a  certain  norm  of  the  solution  is  bounded  by  the  same 
norm  of  the  previous  time  level,  where 

ikiNkl- 

We  notice  that  the  RK2  scheme  (4.5)  above  is  the  same  as  Equation  (3.9),  where 
we  have  already  proven  in  Appendix  C  that  this  time-integrator  is  indeed  second  order 
accurate,  and  from  Chapter  III,  that  the  above  RK2  method,  in  conjunction  with  our 
fourth-order  FD  stencil,  is  consistent  with  the  continuous  PDE  of  our  IVP.  Furthennore, 
we  analyzed  the  stability  of  this  RK2  method  in  the  previous  chapter  using  Von  Neumann 
Stability  Analysis.  Therefore,  if  we  now  substitute  a  linear  combination  of  Euler  steps  as 
defined  by  (4.6),  into  the  above  RK2  method,  we  find  from  [13]: 
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q"'=q-+/(k,+k,)  , 

qm  =qn +\tF{tn,qn)  =  s(At,qn)  , 

=  2^’  +*j)  ’ 

qm  =  q(l)  +  A ,tF(tn,qm)  =  s(At,q{X)) 

1  n  ,  1  (1*) 

=  2*  V  ’ 

n+ 1  1  In,  (1»)\ 

q  =-\q  +q  )  > 

=  ^£(0,<f)  +  ^£(A?,g(1))  , 

where  if  the  forward  Euler  method  is  SSP  under  its  specific  CFL  time-step  restrictions, 
then  Shu  and  Osher  [14]  showed  that  higher-order  methods  constructed  as  linear 
combinations  of  forward  Euler  steps  will  also  be  SSP. 

Now  that  we  have  rewritten  the  RK2  method  as 

qn+'  =^s(0,q")  +  ^£(At,qm)  (4.7) 

the  intent  is  to  use  (4.7)  to  develop  both  a  slow  and  fast  RK2  scheme  to  be  used  within  the 
slow  and  fast  domains  of  our  three  sub-domain  grid  developed  in  the  previous  section. 
We  look  again  to  the  works  by  Constantinescu  and  Sandu  [13],  where  they  show  how  we 
can  extend  the  above  RK2  base  method  to  a  second-order  multi-rate  partitioned  Runge- 
Kutta  (MPRK)  scheme,  where  the  MPRK  scheme  can  be  applied  to  multiple  partitions, 
with  “m”  denoting  the  ratio  between  the  time-steps  associated  with  the  fast  and  slow  sub- 
domains  on  the  same  time  level. 

For  our  analysis,  we  consider  only  the  case  where  our  grid  is  refined  at  most  once, 
such  that  the  grid  consists  of  at  most  two  levels,  which  is  seen  in  Figures  25  -  27. 
Additionally,  for  a  two  level  grid,  one  level  corresponds  to  the  fast  sub-domain,  while  the 
other  level  is  associated  to  the  slow  sub-domain,  such  that  the  fast  domain  is  refined  at  a 
two-to-one  ratio.  Therefore,  we  find  that  if  the  slow  domain  has  elements  of  length  At  , 
than  the  fast  domain  will  have  length  Ax/2.  Furthermore,  the  ratio  between  the  time- 
steps  associated  with  the  fast  and  slow  sub-domains  is  also  two-to-one,  such  that  m  =  2  . 

Fet  us  now  look  at  Figure  28,  where  we  find  the  Butcher  tableau  for  m  =  2  ,  such 
that  the  base  method  as  defined  by  (4.5)  can  be  rewritten  for  the  slow  and  fast  methods, 
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where  the  fast  and  slow  methods  weight  coefficients  are  1  /  mb  repeated  in  times. 
Additionally,  the  slow  method  repeats  the  base  methods  stages  m  times  with  a  time-step 
of  At ,  while  the  fast  method  must  perform  m  steps  of  the  base  method  with  a  time-step  of 
At/m  .  Furthermore,  [13]  shows  how  this  technique  of  partitioning  a  base  RK  method  can 
be  extended  from  m  =  2  to  arbitrary  /«’ s. 
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Figure  28.  Butcher  Tableau  for  RK2  Method  and  its  Associated  Slow  and  Fast  RK2 

Equivalents 


We  have  already  seen  how  we  can  rewrite  the  base  RK2  method  using  a  linear 
combination  of  Euler  steps.  Therefore,  let  us  now  represent  the  slow  and  fast  methods 
shown  in  Table  4  [13]  using  this  same  approach,  where  the  reader  is  encouraged  to  go  to 
[13]  for  proof  of  conservation  and  accuracy  of  these  methods. 

Looking  at  Table  4,  we  notice  that  each  line  in  the  table  is  part  of  an  iterative 
process  for  computing  the  solutions  within  each  sub-domain  of  the  grid.  For  example, 
qF  and  qns  are  the  initial  solutions  on  both  the  slow  and  fast  sub-domains,  which  initially 

correspond  to  the  exact  solution  at  time  tn .  Therefore,  once  we  initialize  the  solution 
vectors  q"F  and  qns ,  we  have  all  the  information  needed  to  move  to  the  second  line  in  the 

table,  and  so  on.  Furthermore,  the  sequence  in  which  each  solution  vector  is  computed, 
resolves  the  issue  for  how  to  compute  the  average  value  for  grid  points  within  the  slow 
domain,  that  are  needed  when  evaluating  our  fourth-order  FD  stencil  in  the  buffer 
regions. 
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Fast  Method 

Slow  Method 

Base  Method 

(buffer  region) 

(slow  region) 

qnF 

q"s 

Ps 

q1f  =£F(f 

Ps  =£s(At,qnF,q'‘) 

qf  =s(At,q'ls) 

qP=sF(%,qV,qf) 

qP=Ss(At,qf,qf) 

qp  =s(At,qf) 

qp=WF+qP) 

qf  =  ql 

qf  =£f  (f>  qf\qns) 

qf  =£S(At,qF),qs) 

qf  =  qf 

q{P=sF(%,q?,q?) 

q(p  =  ss(At,q(Pqp) 

qp  =  qp 

qf  =Ws+WP +WP 

qP=WPqP) 

Table  4.  MPRKo  Algorithm  Used  to  Simultaneously  Solve  the  IVP  on  Both  the  Fast 

and  Slow  Sub-domains 


D.  IMPLEMENTATION,  RESULTS  AND  ANALYSIS 

The  implementation  of  this  iterative  time-integration  process  is  quite  simple.  In 
fact,  the  most  difficult  step  in  advancing  each  slow  and  fast  solution  is  within  the  spatial 
discretization  scheme.  Therefore,  once  the  basic  finite  difference  framework  has  been 
established,  as  was  shown  in  Section  B  of  this  chapter,  the  individual  computations  of 
each  Euler  slow  or  Euler  fast  solution  is  trivial. 

From  Table  4,  we  notice  that  the  base  method  is  equivalent  to  Equation  (4.7),  and 
only  requires  knowledge  of  the  solutions  on  the  purely  slow  domain.  However,  the  fast 
method  performs  four  Euler-fast  steps,  where  each  step  requires  knowledge  from  the  fast 
solution  and  the  slow  solution,  if  solving  for  grid  points  within  the  buffer  region. 
Likewise,  the  slow  method  in  the  buffer  region  also  requires  information  from  the  slow 
solution  and  fast  solution.  In  other  words,  if  a  grid  point  lives  in  either  a  purely  fast  or 
purely  slow  region,  then  the  base  method  will  be  used  for  the  purely  slow  with  only  two 
Euler  slow  steps  each  of  size  At ,  and  the  purely  fast  with  four  Euler  fast  steps  each  of 
size  At/ 2  .  It  is  only  within  the  buffer  region  that  the  MPRK2  method  requires  either  four 
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Euler-fast  or  four  Euler-slow  steps,  each  requiring  information  from  both  the  fast  and 
slow  solutions  to  compute  the  RHS  FD  stencil. 

Let  us  now  look  at  the  results  comparing  the  MPRfG  method  against  the  single¬ 
rate  RK?  for  solving  our  ID,  first  order,  advection  equation  with  constant  wave  speed  and 
periodic  boundary  conditions  for  the  smooth  Gaussian  function  on  the  non-conforming 
three  sub-domain  grid  as  defined  in  Figure  25. 


Figure  29.  RK2  vs.  MPRK2  Convergence  Results  using  Approximately  100  Degrees  of 

Freedom 


From  the  above  figure,  we  notice  that  both  the  MPRK2  and  RK2  results  converge 
on  the  order  of  0(At2)  as  At  — » 0 .  However,  for  step-sizes  At  approximately  less 

than  1 0  3 ,  we  begin  to  notice  the  spatial  error  from  the  fourth-order  FD  stencil. 
Furthermore,  we  also  see  that  for  a  given  time-step  size,  that  the  single-rate  RK2  method 
is  clearly  more  accurate  than  the  MPRK2  for  our  particular  choice  of  spatial  element  sizes 
A*  coarse  =  0.02  and  Ax  fine  =  0.01 .  In  addition,  Figure  30  shows  a  plot  comparing  the 
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efficiency  of  each  method  for  these  same  parameters,  such  that  for  this  size  of  problem, 
i.e.,  total  number  of  grid  points  equal  to  103,  we  notice  that  the  single-rate  Rffi  method  is 
also  more  efficient. 

Given  these  results,  we  now  look  to  increase  the  number  of  points  in  the  spatial 
domain,  to  determine  if/when  the  multi-rate  method  will  be  more  efficient.  Before  we  do 
this,  let  us  first  look  at  the  computational  efficiency  of  each  method,  such  that  the  reason 
for  developing  a  multi-rate  method  in  the  first  place  was  to  have  a  more  efficient 
numerical  code  capable  of  taking  larger  time-steps  within  the  coarse  regions  and  smaller 
time-steps  within  the  fine  regions  of  the  spatial  domain.  As  a  result,  the  multi-rate  method 
should  indeed  be  more  computationally  efficient  than  a  single-rate  code  that  is  restricted 
to  taking  a  time-step  dependant  on  the  smallest  grid  size  element  in  the  spatial  domain. 


Figure  30.  RK2  vs.  MPRK2  Efficiency  Results  using  Approximately  100  Degrees  of 

Freedom 
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Both  the  single-rate  and  multi-rate  methods  are  computed  on  the  same  grid,  such 
that  the  multi-rate  method  solves  for  the  purely  fast  and  purely  slow  solutions  using  the 
same  base  method  as  the  single-rate  method,  except  for  the  purely  slow  solutions  are 
capable  of  taking  a  time-step  twice  as  large  as  the  solution  computed  on  the  fast  domain. 
Therefore,  the  only  difference  in  computations  is  in  the  interface  regions,  which  are 
typically  small  compared  to  the  purely  fast  and  purely  slow  regions.  In  other  words,  we 
truly  have  a  multi-rate  method,  such  that  the  slow  regions  do  in  fact  use  a  step-size  twice 
as  large  as  the  fast  regions. 

Analyzing  the  efficiency  of  our  multi-rate  time-integration  method,  shows  that  the 
speedup  of  this  method  is  equivalent  to  the  ratio  of  the  total  work  done  for  the  single-rate 
scheme  with  its  restrictive  time-step  value  ( At/2  ),  to  the  total  computational  work  of  the 
multi-rate  scheme.  Therefore,  if  we  consider  that  the  multi-rate  method  uses  a  step-size  of 
At/2  on  the  fast  domain,  including  interface  points,  with  NF+N,  grid  points,  and  a 

step-size  of  At  on  the  slow  domain  with  Ns  grid  points,  then  the  multi-rate  speedup  is 


mNr 

m(NF+NI)  +  (m/2)Ns  ’ 


(4.8) 


where 


Nt  =  total  number  of  grid  points  NF  =  number  of  purely  fast  grid  points, 

=  NF  +  Ns  +  Nj,  Nj  =  number  of  interface  grid  points, 

m  =  number  of  fast  steps,  Ns  =  number  of  purely  slow  grid  points. 


Furthennore,  we  find  that  since  the  total  number  of  interface  points  N ,  «;  min(As,  NF) , 
then  for  in  =  4  ,  Equation  (4.8)  simplifies  to 


2  Nt 

2  Nf+Ns 


(4.9) 


Additionally,  taking  the  limit  of  S  as  Ns  — »  NT ,  such  that  AF  — >  0 ,  then  we  find  that 


lim  S  =  2, 

Ns—>Nr 
Nf^>  0 


(4.10) 
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where  Equation  (4.10)  shows  that  the  maximum  theoretical  speedup  of  our  MPRK2 
solution  is  equal  to  2,  regardless  of  the  value  of  m.  This  makes  sense  as  we  have  assumed 
that  there  is  only  one  level  of  refinement  at  a  two-to-one  ratio. 

Below,  Figure  3 1  shows  the  MPRKt  vs.  RK2  wallclock  time  results  versus  time- 
step  size  for  Ax  coarse  =0.02  and  Ax  fine  =0.01.  Using  (4.8),  we  find  that  the 
theoretical  speedup  should  by  approximately  4  percent  for  Nr  =  1 03  grid  points. 

However,  from  this  figure,  we  notice  that  there  was  no  computational  speedup.  In  fact, 
the  single-rate  method  was  still  more  efficient  for  the  above  parameters.  Note,  that  the 
difference  between  the  theoretical  speedup  and  lack  of  computational  speedup  (in  this 
case  specifically)  is  most  likely  attributed  to  other  coding  inefficiencies,  and  not  the 
MPRK2  scheme.  Therefore,  we  expect  to  achieve  close  to  the  theoretical  speedup  as  we 
increase  the  Degrees  of  Freedom  (DoF). 


MPRK2  vs.  RK2  Wallclock  Time  for  dxc  =  0.02  and  dxf  =  0.01 


Figure  3 1 .  RK2  vs.  MPRK2  Effeciency  Results  (-100  DoF) 
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If  we  now  consider  Ax  coarse  =  0.00125  and  Ax  fine  =  0.000625  ,  then  the  total  number 
of  grid  points  is  approximately  1600.  Figure  32  displays  the  convergence  results  for  these 
new  parameters,  such  that  we  find  the  single-rate  RfU  method  is  still  more  accurate  for  a 
given  time-step  size;  however,  Figure  33  shows  that  the  MPRKo  method  is  now  more 
efficient  for  any  given  time-step. 


Figure  32.  RK2  vs.  MPRKo  Convergence  Results  using  Approximately  1600  DoF 


Comparing  the  theoretical  versus  computational  speedup  for  the  newly  defined 
parameters,  Ax  coarse  =0.00125  and  Ax  fine  =0.000625,  then  using  Equation  (4.9), 
we  notice  the  following  results  where  NT  =  1,603 ,  N F  =  5  ,  and  Ns  =  1,596 : 

IN 

S «  T  *  1.996. 

2  Nf  +  Ns 

This  shows  that  for  the  problem  above,  which  has  approximately  1,600  grid  points,  that 
the  theoretical  speedup  should  be  approximately  2  (100%  speedup).  Using  the  data  from 
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Figure  33,  we  find  that  the  computational  speedup  we  achieved  using  the  MPRK2  scheme 
is  approximately  1.67  (67%  speedup).  These  results  are  significant,  and  are  exactly  what 
we  hoped  to  achieve.  In  general,  we  have  proven  that  the  multi-rate  method  formulated  in 
Table  4,  is  indeed  more  efficient  than  the  single-rate  RK2.  Therefore,  given  that  both  the 
RK2  and  MPRK2  time-integration  schemes  are  on  the  order  of  0(At2) ,  the  primary 
consideration  one  must  make  between  the  two  methods  is  whether  one  desires  a  more 
accurate,  but  less  efficient  approximation,  or  a  slightly  less  accurate,  yet  more  efficient 
approximation  to  the  IVP.  In  addition,  let  us  also  look  at  the  estimated  errors  of  each 
method  versus  the  computational  time,  such  that  Figure  34  again  shows  that  the  single¬ 
rate  RK2  method  is  still  more  efficient  overall  in  reaching  a  specific  level  of  accuracy, 
than  the  multi-rate  method. 


MPRK2  vs.  RK2  Wallclock  Time 


Figure  33.  RK2  vs.  MPRK2  Efficiency  Results  using  Approximately  1600  DoF 


Although  Figure  33  demonstrates  that  for  any  given  time-step  the  multi-rate 


method  was  computationally  more  efficient,  the  results  in  Figure  34  were  initially 
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discouraging.  However,  it  can  be  easily  shown  that  the  increased  accuracy  of  the  single¬ 
rate  method  is  attributed  to  the  fixed  non-uniform  grid,  whereas  the  multi-rate  method 
loses  accuracy  due  to  this  same  spatial  mesh.  In  other  words,  the  single-rate  method 
should  indeed  be  more  accurate  as  it  solves  the  IVP  everywhere  on  the  fixed  non-uniform 
grid  with  a  single  time-step  associated  with  the  smallest  element  within  the  spatial 
domain;  whereas,  the  multi-rate  method  uses  a  time-step  twice  as  large  within  the  entire 
slow  sub-domain.  From  the  results  above,  we  see  that  the  slow  sub-domain  consisted  of 
1,596  grid  points  out  of  a  total  1,603  points.  Although  the  initial  grid  was  refined 
underneath  the  bell  of  the  Gaussian  curve,  the  right  moving  wave  was  rarely  inside  the 
fast  region;  therefore,  the  numerical  solution  lost  accuracy  as  the  wave  moved  from  the 
fast  region  into  the  larger  slow  region. 


MPRK2  vs.  RK2  Effeciency  for  dxc  =  0.00125  and  dxf  =  0.000625 
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Figure  34.  RK2  vs.  MPRK2  Efficiency  Results  using  Approximately  1600  DoF 
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V.  SUMMARY  AND  CONCLUSIONS 


Throughout  this  thesis  we  have  shown  how  to  use  the  method  of  finite 
differences,  as  well  as  a  few  explicit  time-integration  schemes,  in  order  to  construct  a 
discretized  fonn  of  our  continuous,  hyperbolic  initial  value  problem,  such  that  we  were 
able  to  successfully  solve  this  IVP  numerically.  In  addition,  we  demonstrated  how  to  use 
a  linear  combination  of  strong  stability  preserving,  explicit  forward  Euler  steps,  in  order 
to  develop  a  2nd  order  Runge-Kutta,  multi-rate  time-integration  scheme,  where  the  results 
from  Chapter  IV  demonstrate  that  this  MPRK2  method  is  indeed  more  computationally 
efficient  than  its  equivalent  single-rate  RKo  method  for  any  given  time-step  size. 

However,  there  is  a  great  deal  of  research  that  must  still  be  accomplished. 
Therefore,  I  recommend  the  following  areas  for  future  research: 

1)  Incorporate  an  adaptive  mesh  refinement  method  in  conjunction  with  the  multi¬ 
rate  scheme. 

2)  Increase  the  dimensionality  of  the  problem. 

3)  Vary  the  initial  and  boundary  conditions. 

4)  Develop  other  multi-rate  time-integration  schemes  using  AB2  or  BDF2. 

Although  the  one-dimensional  efficiency  results  comparing  time-step  sizes  versus  wall- 
clock  time  were  successful,  as  shown  in  Figure  33,  it  would  be  beneficial  to  see  what 
impact  incorporating  recommendation  (1)  would  have  on  improving  the  multi -rate 
accuracy  results  displayed  in  Figure  34.  If  the  refined  sub-domain  within  the  spatial  grid, 
Figure  25,  were  to  adapt  in  time  to  stay  under  the  bell  of  the  Gaussian  curve,  then  we 
should  expect  that  the  MPRK2  method  would  not  only  be  more  efficient  for  a  given  time- 
step  size,  but  also  achieve  a  particular  level  of  accuracy  faster  than  the  single-rate  RK2 
method. 

The  second  recommendation  is  to  increase  the  dimensionality  of  the  problem 
being  solved.  If  we  are  able  to  improve  both  the  efficiency  and  accuracy  of  our  numerical 
solution  in  only  one-dimension,  then  it  follows  that  an  adaptive  multi-rate  time- 
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integration  scheme  should  also  improve  the  speedup  in  higher  dimensions.  In  addition  to 
this  recommendation,  it  would  also  be  beneficial  to  incorporate  other  initial  conditions, 
such  that  we  should  analyze  how  well  the  multi-rate  method  is  capable  of  solving  the 
same  test  case  IVP  with  either  a  square  wave  or  another  Gaussian  wave  with  source. 

Additionally,  future  analysis  should  be  done  on  other  PDE’s  or  systems  of  PDE’s. 
Although  the  1st  order  advection  equation  was  a  good  test  case  PDE,  the  linearized 
shallow  water  equations,  Burgers’  equation,  or  even  Maxwell’s  equations,  each  allow  for 
a  more  thorough  analysis  of  how  well  the  MPRK2  time-integration  scheme  compares  to 
any  other  single-rate  method.  Furthennore,  although  it  is  important  to  analyze  a  particular 
multi-rate  approach  with  its  equivalent  single-rate  method,  as  was  shown  in  this  thesis, 
other  multi-rate  time-integration  schemes  should  also  be  developed  in  order  to  compare 
each  multi-rate  approach  against  another. 
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APPENDIX  A:  TSE  OF  EQUATION  2.1 


We  begin  the  proof  for  Equation  2.1  by  using  Taylor  series  expansion  (TSE)  to 
expand  each  grid  point  about  the  point  qf ,  such  that  for 


dq,  _  g,  -qt- 1 

dx  Ax 


+  0(  Ax) 


We  have  the  following  TSE  for  q" , : 


n  f  K  .ns  n  \  ,  (Ax)2  3  a 

q  .  =  qix,  -  Ax.t  )  =  q - —Ax-\ - + - 0( Ax  ) 

M  '  1  8x  Dx2  2! 


Substituting  q"_x  back  into  Equation  2. 1,  we  then  get 


d_£_ 

dx 


q;~ I  ft'-y-Ar+^^'-OCAr1) 

dx  dx  2 ! 

Ax 


dq"  A  (a*)2  /i/A  3 \ 

Ax  —  ^  +  0(Ax  )  dqn  d2 qn  (Ax) 


dx 


ax2  2! 

Ax 


dx  dx2  2 ! 


+  <9(Ax2) 


0(  Ax) 


dq ; 


dx 


+  0(Ax) 


Therefore,  we  have  shown  that  the  backward  difference  formula  for  the  first  order 
partial  derivative  qx(xi,t" )  is  equal  to  the  exact  solution  plus  an  error  tenn  on  the  order 
of  O(Ax) . 
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APPENDIX  B:  4TH  ORDER,  CENTERED  DIFFERENCE  STENCIL 


We  begin  the  construction  of  this  FD  stencil  by  first  deciding  which  grid  points 
we  want  to  use,  and  then  expanding  about  the  point  (xi,tn)  using  TSE. 

t 


i—  3  i—2  i— 1  i  i+ 1  i+2  i+  3 

Figure  35.  Uniform  Grid  in  Space  and  Time 


From  Figure  35,  we  can  see  that  there  are  many  grid  points  to  choose  from  when 
building  a  new  FD  method.  However,  for  the  construction  of  this  4th  order,  centered 
difference  stencil,  we  will  use  the  following  grid  points: 

*?(•*';+ 1’^ )  ’  ci{xi+2^  )  ’  q (t-i ’ ^ )  ’  q (t-2’ ^ ) 

The  TSE  for  each  point  are 


n  /  .  „  „  dq ”  d2 q”  (Ar)2 

qM=q(xi+Ax,t  )  =  qt  +  —  /Ax  +  —, - — 

ox  ox  21 


d3q"  (Ax)3 


+  0(Ax4) 


n  /  n\  n  dq”  S7tf”(2A>c)^  d3q"  (2/Axf 

qh  =  q(x<  +2Ax,n  =  q”  +~~2Ax+  / 

ox  ox  21  ox  3 ! 


+  0[(  2  Ac)4] 


(B.2) 


77 


q’U  =  q(X[  -Ax,tn)  =  q"-^Ax  +  ^7-  -  ^  ^  +  0( Ax4 ) 


3x 


5x“  2!  ax3  3! 


(5.3) 


ql2  =  q(xi-2Ax,tn)  =  q"-^-2Ax  +  —  4rr  +  0[(2  Ax)4  ] 


dx 


dxz 


2! 


ax3 


3! 


(BA) 


Once  we  have  expanded  each  point,  we  then  need  to  consider  the  following  sum 
of  equations 


A 


+  B 


+  C 


+  D 


Ci=<  + 

n  n  , 

<h+2  =  <1,  + 

M  =  7," 


7,- 2  =  7, 


dx 

n  ,d£ 

dx 

d£_ 
dx 

n  d£ 

dx 


+8X^  +  0(Ax<) 


dx  2! 

^2  n 


dr  3! 


2Ax  +  — -  +  <9[(2Ax)4] 


ar 


2! 


ax3 


3! 


ax2  2!  ax3  3! 


2Ax  +  ^^-^^  +  0[(2Ax)4] 


ax2 


2! 


ax3 


3! 


such  that  we  must  sum  the  columns  of  the  system  to  determine  what  the  coefficients  A,  B, 
C,  and  D  must  be.  We  do  this  by  setting  all  columns  containing  the  derivatives  we  do  not 
want  equal  to  zero.  For  example,  we  are  building  a  FD  stencil  for  the  first  derivative; 
therefore,  the  column  containing  the  first  partial  derivative  will  be  set  equal  to  one  and  all 
others  set  equal  to  zero. 

Summing  the  columns  from  the  equations  above  gives  us: 


Mm  +  Bq'-+ 2  +  C  q"_x  +  Dq"_2 

(5.5) 

q'^A  +  B  +  C  +  D) 

(5.6) 

8c/n 

Hi  Ax(A  +  2B  C  2D) 

dx 

(5.7) 

3  q,2  (A  +  4B  +  C  +  4D) 

dx  2! 

(5.8) 
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8  qt  (Ax)  (A  +  8B_C_8D)  (5.9) 

5x3  3! 

0(Ax4)(A  +  16B  +  C  +  16D)  (5.10) 


Now,  in  order  to  eliminate  tenns,  we  set  equations  (5.5),  (5. 9),  and  ( BAG )  equal 
to  zero  and  Equation  (5. 7)  equal  to  one,  and  then  solve  for  A,  B,  C,  and  D,  such  that: 


(A  +  2B -C -  2D)  =  1 
(A  +  4B  +  C  +  4D)  =  0 
(A  +  8B  -C  -  8D)  =  0 
(A  + 165  +  C  + 165>)  =  0 


125  +  12Z)  =  0  5  =  -£) 

-3A-3C  =  0  ->  4  =  -C 


We  can  then  substitute  these  two  equations  back  into  the  first  four  to  get 


(2A-4D)  =  1  | 
(2A-16D)  =  Oj 


A-2D  =  \\  1  2 

2  \  ->  Z)  =  — ,  ^  =  - 
^4  -  8D  =  0  12  3 


2  1 

Using  these  two  values,  we  then  get  C  =  — ,  and  5  = - . 

3  12 

Another  possible  method  for  quickly  solving  the  above  system  of  equations  is  to 
use  the  method  of  LU  decomposition  by  writing  the  system  in  matrix  form  Ax  =  b ,  and 
then  solving  for  x,  such  that: 


"1 

2 

-1  -2] 

~  A~ 

'll 

A 

1 

4 

1  4 

B 

0 

B 

1 

8 

-1  -8 

C 

0 

9  X  — 

C 

1 

16 

1  16_ 

D 

0_ 

D 

We  have  now  found  the  four  coefficient  values  for  our  4th  order,  centered 
difference  stencil;  therefore,  we  can  now  use  Equations  ( 5.5 )  -  (5.7)  to  construct  the  FD 
stencil. 
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dq’;  (B.5)-(B.6) 

dx  Ax 


dq"  _  C Aq"+l  +  Bq”+2  +  Cq”_x  +  Dq"_2 )  -  (q”  (A  +  B  +  C  +  £>)) 


fix 


Ax 


a?,’ 


1 


1 


T  *7(+l  1  0  ^1+2  a  *7i-l  1  ~  2 

yj  tZ  3  tZ  y 


-(tf(0)) 


fix  Ax 

Finally,  we  have  our  4th  order,  centered  difference  stencil  for  the  first  derivative 


as 


dx  12  Ax 


+  <9(Ax4). 
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APPENDIX  C:  2nd  ORDER  RK2  PROOF 


The  following  is  a  proof,  which  shows  that  regardless  of  the  spatial  discretization 
method,  the  order  of  accuracy  for  any  time-integration  method  is  independent  of  the 
spatial  scheme.  First  let  us  begin  with  Equation  (3.7), 


kx=AtF(tn,q") 

(C.0.1) 

k 2  =  AtF(t',+l ,  q"  +  kx ), 

(C.0.2) 

q"+l  =q"  +^{kl+k2) 

(C.0.3) 

and  the  following  expression, 

q:=F(qn), 

(C.l) 

where  we  can  substitute  C.0.1  and  C.0.2  into  C.0.3  to  arrive  at: 


q"+l=q"  +^-(F{qn)  +  F{q"+MF{q")))  . 

If  we  now  use  Taylor  series  expansion  (TSE)  for  both  the  left  and  right  hand  sides  of  the 
above  expression  for  RK2  we  achieve  the  following: 

A/2  At  f  r)Fn 

q»+q:^  +  q:'+0{AS)  =  q’'+—  F(qn)  +  F(q'‘)  +  ^—AtF(qn)  +  0(At2) 

2  2  y  uc{ 
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q"  =F(q")  +  0(At2), 


(C.3) 


where  we  easily  recognize  that  Equation  (C.3)  is  equivalent  to  (C.l)  with  the  addition  of 
higher  order  terms.  However,  it  must  be  shown  that  the  terms  on  the  order  of  0(At ) 
from  (C.2)  do  not  cancel  each  other  out  when  we  simplified  the  expression  to  Equation 
(C.3). 

Therefore,  if  we  look  at  the  higher  order  terms  on  the  LHS,  we  find  that 

A  p 

0(At3)  =  q;t—  +  0(At4), 

6 

whereas  the  higher  order  tenns  on  the  RHS  are  equivalent  to 

0(At 3)  =  +  0(At4)  . 

4  dq 

If  we  assume  these  two  expressions  are  equal  to  each  other  we  find, 

+  =  +  (C.4) 

6  4  dq 

Now,  if  we  use  the  following  equalities: 


F{qn), 


ay 

dt 2 


dFn 

dq 


F{q\ 


ay 

st 3 


d2Fn  , 

^1- F2(qn)  +  F(q ") 

dq~ 


then  substituting  into  Equation  (C.4), 


F2(qn)  +  F(q ") 


At' 

~6~ 


+  0(At4)± 


At 3  d2Fn 

4  dq2 


F2(q")  +  0(At4), 


we  find  that  these  two  expressions  are  not  equal  to  each  other;  therefore,  showing  that  the 
higher  order  terms  in  Equation  (C.2)  do  not  cancel  each  other,  and  proving  that  by  using 
RK2,  Equation  (C.3)  is  equivalent  to  (C.l)  with  second  order  accuracy. 
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APPENDIX  D:  DERIVATION  OF  AB2  AND  BDF2 


In  order  to  derive  AB2,  we  begin  with 


and  assume  we  want  to  find 


dq" 


dt 


=  F(q‘*i)- 


Using  TSE  to  expand  about  t  2  for  the  solutions  at  q  and  q"  yields 


«’*■ =<r4  +y«r4 +ff--?r4 +o(a»j) 

9 '  =  ?“*  -y  ?r4  +^9,r4  +0(A/5). 


Next,  isolating  g',+2  provides  us  with, 


Atqp  =  q'I+l-q"  +0(  At') 


F 


{«*) 


dq" +2 


dt 


At 


-  +  0(At2)  . 


Next,  let  us  expand  F^q"^  and  F^q"  1  j  about  q"+1  to  arrive  at: 


f  ■  A  t  ^ 


F(q')  =  F  —  ?;*i+0(A/2) 


2 

3  At 


F(qn-l)  =  F  qn+h-— q"^  +  0{At2) 


Now,  taking  a  TSE  about  F^q'" 2  j ,  we  find  that, 
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where  in  order  to  eliminate  O(At)  terms,  requires: 

3  F[qn)-F  [qn-x )  =  IF  (qn+>  j  +  0(At 2 ). 


Therefore,  we  can  finally  show  that, 

F(F-)^F{q-)-l-F(q-')  +  0(M1) 

„'!+1  —  nn  'X  1 

1  +0(At2)  =  ^F(q")-^F(q"-l)  +  0(At2) 

^  =  f  f^-Lf^)^) 

qn+l  =  q"  +  y  (3  F(q’> )  -  F  (qn~l ))  +  0(  At3). 


This  concludes  the  construction  of  ABi  using  TSE  and  matches  Equation  (3.2). 

In  contrast  to  the  AB/>  methods,  which  expand  about  q"+1 ,  the  BDF/>  methods  are 
fully  implicit  with  the  RHS  at  q"+l .  First,  let  us  take 


dqn+1 

dt 


where  if  we  use  TSE  to  expand  q"  and  q"  1  about  qn+1 ,  we  have  the  following: 

At2 

qn  =  qn+1  -  Atqf  +  - +  0(At3) 
qn-]  =  q"+l  -lAtqf  +^^-ql+x  +0{At2). 
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Next,  eliminating  the  terms  O(Ar)  requires  4 q"  -  q"  1 : 


4  qn  -  q'1-'  =  4  q"+l  -  Atq't,+l  +  —  q„  +  0(  AF) 


-  q"+1  -  2Atq"+l  +  q”+l  +  0(Af 3 ) 


4 qn  -qn~l  =  3q  -2Atq';+l  +  0(AF). 


If  we  now  solve  for  q"+1 ,  then 


4q"  -  q"-'  =  3 q"+l  -  2 A tq’t,+l  +  0(AF) 


lnn+1  _?  on 

,n+ 1  _  2  y  ^q  +  2q 


^  = 


-  +  0(Ar). 


Now,  if  we  substitute  the  TSE  for  qn  and  qn  1  into  F  ( q"  j  and  F(t/"  j ,  then  we  have 

*■(«")  =  f  |  -  A +^1«rl  +  0(A<J) 

V  2  7 

F  ( g"'1 )  =  F  qn+l  -  2Atqn;x  +  q[ ;;+1  +  0( At3 ) 


Again,  using  TSE  to  expand  F^q"  j  and  F(g"  1  j  about  F(V'+1) ,  we  then  achieve 

r)Fn+l 

F(q')  =  F(q"‘ )  +  —At qf  +  0( At2 ) 

77’W  +  l 

F  ( <r 1 )  =  F  ( 1 )  +  —  2Atq"+x  +  0(  At2 ), 

where  we  want  to  eliminate  terms  of  O(At) ,  which  requires  2F  (qn  j  -  F(qn  1  j ,  such  that 

F ( q',+1  )  =  2F(q")-F( qn  '  )  +  0(At2). 

Now,  since  we  know  that  F(qn+X^j  =  q"+x ,  we  can  rewrite  the  above  as  follows: 
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+  0(At2 )  =  IF  ( qn )  -  F  ( qnA )  +  0(At2 ) . 


1  nn+l  —  In"  -1-1  n" 

F  ( q"+l )  =  q';+1  =  ^ - 2 ^  2  q 

After  simplifying  the  above  equation,  we  finally  have  the  expression  for  BDF2: 

qn+l  =^qn  ~qn~X  ^{2F(t\qn)-F{tn-\qn~1)), 
which  is  equivalent  to  Equation  (3.6). 

It  is  also  important  to  note  that  neither  derivation  for  AEfi  or  BDF2  required  any 
knowledge  of  how  we  intend  to  solve  for  the  RHS  F(qn],  oxF^qn -1).  In  fact,  we  have 

already  seen  in  Chapter  II,  that  we  can  represent  the  right  hand  side  by  a  wide  variety  of 
difference  schemes,  depending  on  the  order  of  accuracy  we  want  to  achieve  for  the 
spatial  derivative. 
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